!> \brief Collection of routines to compute boundary conditions

MODULE WALL_ROUTINES

USE PRECISION_PARAMETERS
USE GLOBAL_CONSTANTS
USE MESH_POINTERS

IMPLICIT NONE (TYPE,EXTERNAL)
PRIVATE

REAL(EB), POINTER, DIMENSION(:,:) :: PBAR_P
REAL(EB), POINTER, DIMENSION(:,:,:) :: RHOP,HP,UU,VV,WW
REAL(EB), POINTER, DIMENSION(:,:,:,:) :: ZZP

PUBLIC WALL_BC,TGA_ANALYSIS,HT3D_TEMPERATURE_EXCHANGE

CONTAINS


!> \brief Main control routine for applying boundary conditions.
!>
!> \param T Current time (s)
!> \param DT Current time step (s)
!> \param NM Mesh number

SUBROUTINE WALL_BC(T,DT,NM)

USE COMP_FUNCTIONS, ONLY: CURRENT_TIME
USE CC_SCALARS, ONLY: CFACE_THERMAL_GASVARS
USE TURBULENCE, ONLY: WALL_MODEL
REAL(EB) :: TNOW
REAL(EB), INTENT(IN) :: T,DT
INTEGER, INTENT(IN) :: NM
LOGICAL :: CALL_HT_1D
REAL(EB) :: DT_BC,SLIP_COEF
INTEGER :: IW,IP,ICF,ITW
TYPE(WALL_TYPE), POINTER :: WC
TYPE(SURFACE_TYPE), POINTER :: SF
TYPE(CFACE_TYPE), POINTER :: CFA
TYPE(LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP
TYPE(LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC
TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC
TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1
TYPE(BOUNDARY_PROP2_TYPE), POINTER :: B2

IF (LEVEL_SET_MODE==1) RETURN  ! No need for boundary conditions if the simulation is uncoupled fire spread only

TNOW=CURRENT_TIME()

CALL POINT_TO_MESH(NM)

IF (PREDICTOR) THEN
   UU => US
   VV => VS
   WW => WS
   RHOP => RHOS
   ZZP  => ZZS
   PBAR_P => PBAR_S
   HP     => HS
ELSE
   UU => U
   VV => V
   WW => W
   RHOP => RHO
   ZZP  => ZZ
   PBAR_P => PBAR
   HP     => H
ENDIF

! For thermally-thick boundary conditions, set the flag, CALL_HT_1D, to call the subroutine SOLID_HEAT_TRANSFER.
! This routine is called every WALL_INCREMENT time steps. Also, if any 3D heat transfer calculations are done, set a new "sweep
! direction" (1, 2, or 3)

CALL_HT_1D = .FALSE.

IF (.NOT.INITIALIZATION_PHASE .AND. CORRECTOR) THEN
   IF (WALL_COUNTER==WALL_INCREMENT) THEN
      DT_BC    = T - BC_CLOCK
      BC_CLOCK = T
      CALL_HT_1D = .TRUE.
      HT_3D_SWEEP_DIRECTION = HT_3D_SWEEP_DIRECTION + 1
      IF (HT_3D_SWEEP_DIRECTION>3) HT_3D_SWEEP_DIRECTION = 1
   ENDIF
ENDIF

!$OMP PARALLEL PRIVATE(IW,WC,SF,BC,B1,B2,LP,LPC,CFA,IP,ICF,SLIP_COEF)

! Sweep through all WALL cells and apply boundary conditions

!$OMP DO SCHEDULE(DYNAMIC)
WALL_CELL_LOOP: DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS

   WC=>WALL(IW)
   BC => BOUNDARY_COORD(WC%BC_INDEX)
   B1 => BOUNDARY_PROP1(WC%B1_INDEX)
   IF (IW<=N_EXTERNAL_WALL_CELLS) CALL ASSIGN_GHOST_VALUE(IW,BC,B1)
   IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE WALL_CELL_LOOP
   B2 => BOUNDARY_PROP2(WC%B2_INDEX)
   SF => SURFACE(WC%SURF_INDEX)

   CALL NEAR_SURFACE_GAS_VARIABLES(T,SF,BC,B1,WALL_INDEX=IW)

   IF (.NOT.SF%THERMAL_BC_INDEX==THERMALLY_THICK) THEN
      CALL SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX=IW)
   ELSEIF (CALL_HT_1D) THEN
      CALL SOLID_HEAT_TRANSFER(NM,T,SF%HT_DIM*DT_BC,WALL_INDEX=IW)
   ENDIF

   IF (N_TRACKED_SPECIES>1 .AND. WC%BOUNDARY_TYPE/=OPEN_BOUNDARY .AND. WC%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY) THEN
      CALL CALCULATE_RHO_D_F(B1,BC,WALL_INDEX=IW)
   ENDIF

   IF (WC%BOUNDARY_TYPE==SOLID_BOUNDARY .AND. &
       (ANY(SPECIES_MIXTURE%CONDENSATION_SMIX_INDEX>0) .OR. DEPOSITION .OR. OUTPUT_WALL_QUANTITIES)) THEN
      CALL WALL_MODEL(SLIP_COEF,B2%U_TAU,B2%Y_PLUS,MU_DNS(BC%IIG,BC%JJG,BC%KKG)/RHO(BC%IIG,BC%JJG,BC%KKG),SF%ROUGHNESS,&
                      0.5_EB/B1%RDN,B1%U_TANG)
   ENDIF

   IF (DEPOSITION .AND. .NOT.INITIALIZATION_PHASE .AND. CORRECTOR .AND. .NOT.SOLID_PHASE_ONLY) THEN
      IF (WC%BOUNDARY_TYPE==SOLID_BOUNDARY .AND. B1%NODE_INDEX==0 .AND. ABS(SF%VEL)<TWO_EPSILON_EB .AND. &
          ANY(ABS(SF%MASS_FLUX)<TWO_EPSILON_EB) .AND. ABS(SF%VOLUME_FLOW)<TWO_EPSILON_EB) THEN
         CALL CALC_DEPOSITION(DT,BC,B1,B2,WALL_INDEX=IW)
      ENDIF
   ENDIF

   IF (HVAC_SOLVE .AND. .NOT.INITIALIZATION_PHASE) THEN
      IF (B1%NODE_INDEX/=0) CALL CALC_HVAC_BC(BC,B1,SF)
   ENDIF

   IF (WC%BOUNDARY_TYPE/=OPEN_BOUNDARY .AND. WC%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY) THEN
      CALL CALCULATE_ZZ_F(T,DT,WALL_INDEX=IW)
   ENDIF

   IF (WC%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY) CALL CALCULATE_RHO_F(BC,B1,WALL_INDEX=IW)

ENDDO WALL_CELL_LOOP
!$OMP END DO

! Do lateral heat transfer in thin obstructions

IF (CALL_HT_1D) THEN
   !$OMP DO SCHEDULE(DYNAMIC)
   DO ITW=1,N_THIN_WALL_CELLS
      CALL SOLID_HEAT_TRANSFER(NM,T,3._EB*DT_BC,THIN_WALL_INDEX=ITW)
   ENDDO
   !$OMP END DO
ENDIF

! Loop through all CFACEs and apply heat transfer BCs

!$OMP DO SCHEDULE(DYNAMIC)
CFACE_LOOP: DO ICF=INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS
   CFA => CFACE(ICF)
   IF (CFA%BOUNDARY_TYPE==NULL_BOUNDARY) CYCLE CFACE_LOOP
   SF => SURFACE(CFA%SURF_INDEX)
   BC => BOUNDARY_COORD(CFA%BC_INDEX)
   B1 => BOUNDARY_PROP1(CFA%B1_INDEX)
   B2 => BOUNDARY_PROP2(CFA%B2_INDEX)
   CALL CFACE_THERMAL_GASVARS(ICF,SF,B1,T)
   IF (.NOT.SF%THERMAL_BC_INDEX==THERMALLY_THICK) THEN
      CALL SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,CFACE_INDEX=ICF)
   ELSEIF (CALL_HT_1D) THEN
      CALL SOLID_HEAT_TRANSFER(NM,T,DT_BC,CFACE_INDEX=ICF)
   ENDIF

   IF (N_TRACKED_SPECIES>1) CALL CALCULATE_RHO_D_F(B1,BC,CFACE_INDEX=ICF)

   IF (CFA%BOUNDARY_TYPE==SOLID_BOUNDARY .AND. &
      (ANY(SPECIES_MIXTURE%CONDENSATION_SMIX_INDEX>0) .OR. DEPOSITION .OR. OUTPUT_WALL_QUANTITIES)) THEN
      CALL WALL_MODEL(SLIP_COEF,B2%U_TAU,B2%Y_PLUS,MU_DNS(BC%IIG,BC%JJG,BC%KKG)/RHO(BC%IIG,BC%JJG,BC%KKG),SF%ROUGHNESS,&
                      0.5_EB/B1%RDN,B1%U_TANG)
   ENDIF

   IF (DEPOSITION .AND. .NOT.INITIALIZATION_PHASE .AND. CORRECTOR .AND. .NOT.SOLID_PHASE_ONLY) THEN
      IF (CFA%BOUNDARY_TYPE==SOLID_BOUNDARY .AND. B1%NODE_INDEX==0 .AND. ABS(SF%VEL)<TWO_EPSILON_EB .AND. &
          ANY(ABS(SF%MASS_FLUX)<TWO_EPSILON_EB) .AND. ABS(SF%VOLUME_FLOW)<TWO_EPSILON_EB) THEN
         CALL CALC_DEPOSITION(DT,BC,B1,B2,CFACE_INDEX=ICF)
      ENDIF
   ENDIF

   IF (HVAC_SOLVE .AND. .NOT.INITIALIZATION_PHASE) THEN
      IF (B1%NODE_INDEX/=0) CALL CALC_HVAC_BC(BC,B1,SF)
   ENDIF

   CALL CALCULATE_ZZ_F(T,DT,CFACE_INDEX=ICF)

   CALL CALCULATE_RHO_F(BC,B1,CFACE_INDEX=ICF)

ENDDO CFACE_LOOP
!$OMP END DO

! Apply boundary conditions to SOLID PARTICLES

IF (SOLID_PARTICLES) THEN

   !$OMP DO SCHEDULE(DYNAMIC)
   PARTICLE_LOOP: DO IP=1, NLP

      LP => LAGRANGIAN_PARTICLE(IP)
      LPC => LAGRANGIAN_PARTICLE_CLASS(LP%CLASS_INDEX)
      IF (.NOT.LPC%SOLID_PARTICLE .AND. .NOT.LPC%MASSLESS_TARGET) CYCLE PARTICLE_LOOP
      SF => SURFACE(LPC%SURF_INDEX)
      BC => BOUNDARY_COORD(LP%BC_INDEX)
      B1 => BOUNDARY_PROP1(LP%B1_INDEX)

      CALL NEAR_SURFACE_GAS_VARIABLES(T,SF,BC,B1,LP=LP,PARTICLE_INDEX=IP)

      IF (.NOT.SF%THERMAL_BC_INDEX==THERMALLY_THICK) THEN
         CALL SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,PARTICLE_INDEX=IP)
      ELSEIF (CALL_HT_1D) THEN
         CALL SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX=IP)
      ENDIF

      IF (LPC%SOLID_PARTICLE) THEN
         CALL CALCULATE_ZZ_F(T,DT,PARTICLE_INDEX=IP)
         IF (CORRECTOR) CALL DEPOSIT_PARTICLE_MASS(NM,LP,LPC)  ! Add the particle off-gas to the gas phase mesh
      ENDIF

   ENDDO PARTICLE_LOOP
   !$OMP END DO

ENDIF

!$OMP END PARALLEL

T_USED(6)=T_USED(6)+CURRENT_TIME()-TNOW
END SUBROUTINE WALL_BC


SUBROUTINE ASSIGN_GHOST_VALUE(WALL_INDEX,BC,B1)

USE PHYSICAL_FUNCTIONS, ONLY : GET_SPECIFIC_GAS_CONSTANT
INTEGER, INTENT(IN) :: WALL_INDEX
REAL(EB) :: ARO,ZZ_GET(1:N_TRACKED_SPECIES),RHO_OTHER,RHO_ZZ_OTHER(1:N_TOTAL_SCALARS)
INTEGER :: IIO,JJO,KKO
REAL(EB), POINTER, DIMENSION(:,:,:,:) :: OM_ZZP
REAL(EB), POINTER, DIMENSION(:,:,:) :: OM_RHOP
TYPE(EXTERNAL_WALL_TYPE), POINTER :: EWC
TYPE(OMESH_TYPE), POINTER :: OM
TYPE(MESH_TYPE), POINTER :: MM
TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC
TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1

EWC => EXTERNAL_WALL(WALL_INDEX)
IF (EWC%NOM==0) RETURN

OM => OMESH(EWC%NOM)
IF (PREDICTOR) THEN
   OM_RHOP => OM%RHOS
   OM_ZZP => OM%ZZS
ELSE
   OM_RHOP => OM%RHO
   OM_ZZP => OM%ZZ
ENDIF
MM => MESHES(EWC%NOM)

! Compute average values of RHO and RHO_ZZ at mesh boundary

RHO_OTHER    = 0._EB
RHO_ZZ_OTHER = 0._EB

DO KKO=EWC%KKO_MIN,EWC%KKO_MAX
   DO JJO=EWC%JJO_MIN,EWC%JJO_MAX
      DO IIO=EWC%IIO_MIN,EWC%IIO_MAX
         SELECT CASE(BC%IOR)
            CASE( 1) ; ARO = MIN(1._EB , (MM%DY(JJO)*MM%DZ(KKO))/(DY(BC%JJ)*DZ(BC%KK)) )
            CASE(-1) ; ARO = MIN(1._EB , (MM%DY(JJO)*MM%DZ(KKO))/(DY(BC%JJ)*DZ(BC%KK)) )
            CASE( 2) ; ARO = MIN(1._EB , (MM%DX(IIO)*MM%DZ(KKO))/(DX(BC%II)*DZ(BC%KK)) )
            CASE(-2) ; ARO = MIN(1._EB , (MM%DX(IIO)*MM%DZ(KKO))/(DX(BC%II)*DZ(BC%KK)) )
            CASE( 3) ; ARO = MIN(1._EB , (MM%DX(IIO)*MM%DY(JJO))/(DX(BC%II)*DY(BC%JJ)) )
            CASE(-3) ; ARO = MIN(1._EB , (MM%DX(IIO)*MM%DY(JJO))/(DX(BC%II)*DY(BC%JJ)) )
         END SELECT
         RHO_OTHER = RHO_OTHER + ARO*OM_RHOP(IIO,JJO,KKO)
         RHO_ZZ_OTHER(1:N_TOTAL_SCALARS) = RHO_ZZ_OTHER(1:N_TOTAL_SCALARS) + &
                                           ARO*OM_RHOP(IIO,JJO,KKO)*OM_ZZP(IIO,JJO,KKO,1:N_TOTAL_SCALARS)
      ENDDO
   ENDDO
ENDDO

! Assign density (RHOP), mass fraction (ZZP), and temperature (TMP) to ghost cells at exterior of the mesh

RHOP(BC%II,BC%JJ,BC%KK) = RHO_OTHER
ZZP(BC%II,BC%JJ,BC%KK,1:N_TOTAL_SCALARS) = MAX(0._EB,MIN(1._EB,RHO_ZZ_OTHER(1:N_TOTAL_SCALARS)/RHO_OTHER))
ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(BC%II,BC%JJ,BC%KK,1:N_TRACKED_SPECIES)
CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET,RSUM(BC%II,BC%JJ,BC%KK))
TMP(BC%II,BC%JJ,BC%KK) = PBAR_P(BC%KK,B1%PRESSURE_ZONE)/(RSUM(BC%II,BC%JJ,BC%KK)*RHOP(BC%II,BC%JJ,BC%KK))

END SUBROUTINE ASSIGN_GHOST_VALUE


!> \brief Set various near-surface variables to be used for boundary conditions
!> \param T Time (s)
!> \param SF Pointer to SURFACE derived type
!> \param BC Pointer to BOUNDARY_COORD derived type
!> \param B1 Pointer to BOUNDARY_PROP1 derived type

SUBROUTINE NEAR_SURFACE_GAS_VARIABLES(T,SF,BC,B1,LP,WALL_INDEX,PARTICLE_INDEX)

USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP
REAL(EB), INTENT(IN) :: T
INTEGER, INTENT(IN), OPTIONAL :: WALL_INDEX,PARTICLE_INDEX
TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1
TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC
TYPE(SURFACE_TYPE), POINTER :: SF
TYPE(LAGRANGIAN_PARTICLE_TYPE), POINTER, OPTIONAL :: LP
REAL(EB) :: TSI,RAMP_FACTOR,UBAR,VBAR,WBAR

IF (PRESENT(WALL_INDEX)) THEN
   IF (ABS(SF%T_IGN-T_BEGIN)<=SPACING(SF%T_IGN) .AND. SF%RAMP(TIME_VELO)%INDEX>=1) THEN
      TSI = T
   ELSE
      TSI = T-SF%T_IGN
   ENDIF
   RAMP_FACTOR = EVALUATE_RAMP(TSI,SF%RAMP(TIME_VELO)%INDEX,TAU=SF%RAMP(TIME_VELO)%TAU)
   SELECT CASE(BC%IOR)
      CASE(1,-1)
         VBAR = 0.5_EB*(VV(BC%IIG,BC%JJG,BC%KKG)+VV(BC%IIG,BC%JJG-1,BC%KKG)) - SF%VEL_T(1)*RAMP_FACTOR
         WBAR = 0.5_EB*(WW(BC%IIG,BC%JJG,BC%KKG)+WW(BC%IIG,BC%JJG,BC%KKG-1)) - SF%VEL_T(2)*RAMP_FACTOR
         B1%U_TANG = SQRT(VBAR**2+WBAR**2)
      CASE(2,-2)
         UBAR = 0.5_EB*(UU(BC%IIG,BC%JJG,BC%KKG)+UU(BC%IIG-1,BC%JJG,BC%KKG)) - SF%VEL_T(1)*RAMP_FACTOR
         WBAR = 0.5_EB*(WW(BC%IIG,BC%JJG,BC%KKG)+WW(BC%IIG,BC%JJG,BC%KKG-1)) - SF%VEL_T(2)*RAMP_FACTOR
         B1%U_TANG = SQRT(UBAR**2+WBAR**2)
      CASE(3,-3)
         UBAR = 0.5_EB*(UU(BC%IIG,BC%JJG,BC%KKG)+UU(BC%IIG-1,BC%JJG,BC%KKG)) - SF%VEL_T(1)*RAMP_FACTOR
         VBAR = 0.5_EB*(VV(BC%IIG,BC%JJG,BC%KKG)+VV(BC%IIG,BC%JJG-1,BC%KKG)) - SF%VEL_T(2)*RAMP_FACTOR
         B1%U_TANG = SQRT(UBAR**2+VBAR**2)
   END SELECT
   B1%U_IMPACT = SQRT(2._EB*MAX(HP(BC%IIG,BC%JJG,BC%KKG),TWO_EPSILON_EB))
ELSEIF (PRESENT(PARTICLE_INDEX)) THEN
   UBAR = 0.5_EB*(UU(BC%IIG,BC%JJG,BC%KKG)+UU(BC%IIG-1,BC%JJG,BC%KKG)) - LP%U
   VBAR = 0.5_EB*(VV(BC%IIG,BC%JJG,BC%KKG)+VV(BC%IIG,BC%JJG-1,BC%KKG)) - LP%V
   WBAR = 0.5_EB*(WW(BC%IIG,BC%JJG,BC%KKG)+WW(BC%IIG,BC%JJG,BC%KKG-1)) - LP%W
   B1%U_TANG = SQRT(UBAR**2+VBAR**2+WBAR**2)
ENDIF

IF (SF%TMP_GAS_FRONT > 0._EB) THEN
   B1%TMP_G = TMPA + EVALUATE_RAMP(T-T_BEGIN,SF%RAMP(TIME_TGF)%INDEX)*(SF%TMP_GAS_FRONT-TMPA)
ELSE
   B1%TMP_G =  TMP(BC%IIG,BC%JJG,BC%KKG)
ENDIF
B1%RHO_G = RHOP(BC%IIG,BC%JJG,BC%KKG)
B1%ZZ_G(1:N_TRACKED_SPECIES) = ZZP(BC%IIG,BC%JJG,BC%KKG,1:N_TRACKED_SPECIES)

END SUBROUTINE NEAR_SURFACE_GAS_VARIABLES


!> \brief Calculate the surface temperature TMP_F
!> \param NM Mesh number
!> \param T Time (s)
!> \param WALL_INDEX Optional WALL cell index
!> \param CFACE_INDEX Optional immersed boundary (CFACE) index
!> \param PARTICLE_INDEX Optional Lagrangian particle index
!> \details Calculate the surface temperature TMP_F of either a WALL cell, immersed CFACE cell, or a Lagrangian particle.

SUBROUTINE SURFACE_HEAT_TRANSFER(NM,T,SF,BC,B1,WALL_INDEX,CFACE_INDEX,PARTICLE_INDEX)

USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP,INTERPOLATE1D_UNIFORM, GET_SCALAR_FACE_VALUE
USE COMPLEX_GEOMETRY, ONLY : CC_CGSC, CC_SOLID
USE PHYSICAL_FUNCTIONS, ONLY : GET_SPECIFIC_GAS_CONSTANT,GET_VISCOSITY
USE DEVICE_VARIABLES, ONLY : PROPERTY,PROPERTY_TYPE

REAL(EB), INTENT(IN) :: T
INTEGER, INTENT(IN) :: NM
INTEGER, INTENT(IN), OPTIONAL :: WALL_INDEX,CFACE_INDEX,PARTICLE_INDEX
REAL(EB) :: ARO,QNET,RAMP_FACTOR,RHO_G_2,RSUM_F,PBAR_F,TSI,UN, &
            RHO_ZZ_F(1:N_TOTAL_SCALARS),ZZ_GET(1:N_TRACKED_SPECIES), &
            RHO_OTHER,RHO_OTHER_2,RHO_ZZ_OTHER(1:N_TOTAL_SCALARS),RHO_ZZ_OTHER_2,RHO_ZZ_G,RHO_ZZ_G_2, &
            DDO,PBAR_G,PBAR_OTHER,DENOM,D_Z_N(0:I_MAX_TEMP),D_Z_G,D_Z_OTHER,TMP_OTHER, &
            MU_DNS_G,MU_DNS_OTHER,MU_OTHER,RHO_D,RHO_D_TURB,RHO_D_DZDN,RHO_D_DZDN_OTHER,RSUM_OTHER,SF_HTC,&
            BBB,CCC,PPP,QQQ,RRR,UUU,YYY,WWW,HTC_OLD,RSC_LOC,DTMP,RSUM_G,MU_G
LOGICAL :: SECOND_ORDER_INTERPOLATED_BOUNDARY,ATMOSPHERIC_INTERPOLATION,CC_SOLID_FLAG
INTEGER :: IIO,JJO,KKO,N,ADCOUNT,ICG,ICO
REAL(EB), POINTER, DIMENSION(:,:,:,:) :: OM_ZZP
REAL(EB), POINTER, DIMENSION(:,:,:) :: OM_MUP
REAL(EB), DIMENSION(0:3,0:3,0:3) :: U_TEMP,Z_TEMP,F_TEMP
TYPE(EXTERNAL_WALL_TYPE), POINTER :: EWC
TYPE(OMESH_TYPE), POINTER :: OM
TYPE(MESH_TYPE), POINTER :: MM
TYPE(VENTS_TYPE), POINTER :: VT
TYPE(WALL_TYPE), POINTER :: WC
TYPE(CFACE_TYPE), POINTER :: CFA
TYPE(SURFACE_TYPE), POINTER :: SF
TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1
TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC
TYPE(LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP
TYPE(LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC
REAL(EB), POINTER, DIMENSION(:,:,:) :: OM_RHOP
TYPE(PROPERTY_TYPE), POINTER :: PY

IF (PRESENT(WALL_INDEX)) THEN
   RSUM_G  = RSUM(BC%IIG,BC%JJG,BC%KKG)
   MU_G    = MU(BC%IIG,BC%JJG,BC%KKG)
ELSEIF (PRESENT(CFACE_INDEX)) THEN
   CFA=> CFACE(CFACE_INDEX)
   RSUM_G  = CFA%RSUM_G
   MU_G    = CFA%MU_G
ELSEIF (PRESENT(PARTICLE_INDEX)) THEN
   LP => LAGRANGIAN_PARTICLE(PARTICLE_INDEX)
   LPC => LAGRANGIAN_PARTICLE_CLASS(LP%CLASS_INDEX)
   RSUM_G = RSUM(BC%IIG,BC%JJG,BC%KKG)
   MU_G   = MU(BC%IIG,BC%JJG,BC%KKG)
   IF (LPC%MASSLESS_TARGET) THEN
      PY => PROPERTY(LP%PROP_INDEX)
      IF (PY%HEAT_TRANSFER_COEFFICIENT>0._EB) THEN
         B1%HEAT_TRANS_COEF = PY%HEAT_TRANSFER_COEFFICIENT
      ELSE
         IF (SF%H_FIXED >= 0._EB) THEN
            SF_HTC = SF%H_FIXED
            IF (SF%RAMP_H_FIXED_INDEX > 0) SF_HTC = SF_HTC * EVALUATE_RAMP(T-T_BEGIN,SF%RAMP_H_FIXED_INDEX)
         ELSE
            SF_HTC = -1._EB
         ENDIF
         B1%HEAT_TRANS_COEF = HEAT_TRANSFER_COEFFICIENT(NM,B1%TMP_G-PY%GAUGE_TEMPERATURE,SF_HTC,SF,&
                              PARTICLE_INDEX_IN=PARTICLE_INDEX)
      ENDIF
      B1%Q_CON_F = B1%HEAT_TRANS_COEF*(B1%TMP_G-PY%GAUGE_TEMPERATURE)
      RETURN
   ENDIF
ENDIF

! Compute surface temperature, TMP_F, and convective heat flux, Q_CON_F, for various boundary conditions

METHOD_OF_HEAT_TRANSFER: SELECT CASE(SF%THERMAL_BC_INDEX)

   CASE (NO_CONVECTION) METHOD_OF_HEAT_TRANSFER

      B1%TMP_F  = B1%TMP_G
      IF (RADIATION) B1%Q_RAD_OUT = SIGMA*B1%EMISSIVITY*B1%TMP_F**4

   CASE (INFLOW_OUTFLOW) METHOD_OF_HEAT_TRANSFER  ! Only for WALL cells

      ! Base inflow/outflow decision on velocity component with same predictor/corrector attribute

      SELECT CASE(BC%IOR)
         CASE( 1); UN =  UU(BC%II  ,BC%JJ  ,BC%KK  )
         CASE(-1); UN = -UU(BC%II-1,BC%JJ  ,BC%KK  )
         CASE( 2); UN =  VV(BC%II  ,BC%JJ  ,BC%KK  )
         CASE(-2); UN = -VV(BC%II  ,BC%JJ-1,BC%KK  )
         CASE( 3); UN =  WW(BC%II  ,BC%JJ  ,BC%KK  )
         CASE(-3); UN = -WW(BC%II  ,BC%JJ  ,BC%KK-1)
         CASE DEFAULT; UN = 0._EB
      END SELECT

      IF (UN>TWO_EPSILON_EB) THEN  ! Assume the flow is coming into the domain
         B1%TMP_F = TMP_0(BC%KK)
         WC => WALL(WALL_INDEX)
         IF (WC%VENT_INDEX>0) THEN
            VT => VENTS(WC%VENT_INDEX)
            IF (VT%TMP_EXTERIOR>0._EB) THEN
               TSI = T - T_BEGIN
               B1%TMP_F = TMP_0(BC%KK) + EVALUATE_RAMP(TSI,VT%TMP_EXTERIOR_RAMP_INDEX)*(VT%TMP_EXTERIOR-TMP_0(BC%KK))
            ENDIF
         ENDIF
         B1%ZZ_F(1:N_TRACKED_SPECIES)=SPECIES_MIXTURE(1:N_TRACKED_SPECIES)%ZZ0
      ELSE
         B1%TMP_F = B1%TMP_G
         B1%ZZ_F(1:N_TRACKED_SPECIES) = B1%ZZ_G(1:N_TRACKED_SPECIES)
      ENDIF

      ! Avoid large fluxes at open downwind boundaries

      IF (OPEN_WIND_BOUNDARY .AND. DOT_PRODUCT(BC%NVEC,(/U_WIND(BC%KK),V_WIND(BC%KK),W_WIND(BC%KK)/))<-TWO_EPSILON_EB) THEN
         B1%TMP_F = B1%TMP_G
         B1%ZZ_F(1:N_TRACKED_SPECIES) = B1%ZZ_G(1:N_TRACKED_SPECIES)
      ENDIF

      ! Ghost cell values

      TMP(BC%II,BC%JJ,BC%KK) = B1%TMP_F
      ZZP(BC%II,BC%JJ,BC%KK,1:N_TRACKED_SPECIES) = B1%ZZ_F(1:N_TRACKED_SPECIES)
      ZZ_GET(1:N_TRACKED_SPECIES) = MAX(0._EB,ZZP(BC%II,BC%JJ,BC%KK,1:N_TRACKED_SPECIES))
      CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET,RSUM(BC%II,BC%JJ,BC%KK))
      RHOP(BC%II,BC%JJ,BC%KK) = PBAR_P(BC%KK,B1%PRESSURE_ZONE)/(RSUM(BC%II,BC%JJ,BC%KK)*TMP(BC%II,BC%JJ,BC%KK))

      B1%Q_CON_F = 2._EB*B1%K_G*(B1%TMP_G-B1%TMP_F)*B1%RDN
      IF (RADIATION) B1%Q_RAD_OUT = SIGMA*B1%EMISSIVITY*B1%TMP_F**4

   CASE (SPECIFIED_TEMPERATURE) METHOD_OF_HEAT_TRANSFER

      IF (ABS(B1%T_IGN-T_BEGIN) <= SPACING(B1%T_IGN) .AND. SF%RAMP(TIME_TEMP)%INDEX>=1) THEN
         TSI = T
      ELSE
         TSI = T - B1%T_IGN
      ENDIF

      IF (B1%U_NORMAL>TWO_EPSILON_EB) THEN
         B1%TMP_F = B1%TMP_G
      ELSEIF (SF%TMP_FRONT>0._EB) THEN
         B1%TMP_F = TMP_0(BC%KKG) + &
                    EVALUATE_RAMP(TSI,SF%RAMP(TIME_TEMP)%INDEX,TAU=SF%RAMP(TIME_TEMP)%TAU)*(SF%TMP_FRONT-TMP_0(BC%KKG))
      ELSE
         B1%TMP_F = TMP_0(BC%KKG)
      ENDIF

      DTMP = B1%TMP_G - B1%TMP_F

      IF (SF%H_FIXED >= 0._EB) THEN
         SF_HTC = SF%H_FIXED
         IF (SF%RAMP_H_FIXED_INDEX > 0) SF_HTC = SF_HTC * EVALUATE_RAMP(T-T_BEGIN,SF%RAMP_H_FIXED_INDEX)
      ELSE
         SF_HTC = -1._EB
      ENDIF

      IF (PRESENT(WALL_INDEX)) THEN
         B1%HEAT_TRANS_COEF = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC,SF,WALL_INDEX_IN=WALL_INDEX)
      ELSEIF (PRESENT(CFACE_INDEX)) THEN
         B1%HEAT_TRANS_COEF = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC,SF,CFACE_INDEX_IN=CFACE_INDEX)
      ELSEIF (PRESENT(PARTICLE_INDEX)) THEN
         B1%HEAT_TRANS_COEF = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC,SF,PARTICLE_INDEX_IN=PARTICLE_INDEX)
      ENDIF
      B1%Q_CON_F = B1%HEAT_TRANS_COEF*DTMP
      IF (RADIATION) B1%Q_RAD_OUT = SIGMA*B1%EMISSIVITY*B1%TMP_F**4

   CASE (CONVECTIVE_FLUX_BC,NET_FLUX_BC) METHOD_OF_HEAT_TRANSFER

      IF (ABS(B1%T_IGN-T_BEGIN)<= SPACING(B1%T_IGN) .AND. SF%RAMP(TIME_HEAT)%INDEX>=1) THEN
         TSI = T
      ELSE
         TSI = T - B1%T_IGN
      ENDIF
      HTC_OLD = B1%HEAT_TRANS_COEF
      RAMP_FACTOR = EVALUATE_RAMP(TSI,SF%RAMP(TIME_HEAT)%INDEX,TAU=SF%RAMP(TIME_HEAT)%TAU)
      IF (SF%SET_H) THEN
         B1%Q_CON_F = -RAMP_FACTOR*SF%CONVECTIVE_HEAT_FLUX*B1%AREA_ADJUST
         B1%HEAT_TRANS_COEF = B1%Q_CON_F/(B1%TMP_G-B1%TMP_F+TWO_EPSILON_EB)
      ELSE
         IF (SF%THERMAL_BC_INDEX==NET_FLUX_BC) THEN
            QNET = -RAMP_FACTOR*SF%NET_HEAT_FLUX*B1%AREA_ADJUST
         ELSE
            QNET = -RAMP_FACTOR*SF%CONVECTIVE_HEAT_FLUX*B1%AREA_ADJUST
         ENDIF

         IF (SF%H_FIXED >= 0._EB) THEN
            SF_HTC = SF%H_FIXED
            IF (SF%RAMP_H_FIXED_INDEX > 0) SF_HTC = SF_HTC * EVALUATE_RAMP(T-T_BEGIN,SF%RAMP_H_FIXED_INDEX)
         ELSE
            SF_HTC = -1._EB
         ENDIF
         IF (ABS(SF_HTC) < TWO_EPSILON_EB) THEN
            IF (RADIATION) B1%TMP_F = ((-QNET + B1%Q_RAD_IN)/(B1%EMISSIVITY * SIGMA))**0.25_EB
            B1%HEAT_TRANS_COEF = 0._EB
            B1%Q_CON_F = 0._EB
         ELSE
            ADCOUNT = 0
            ADLOOP: DO
               ADCOUNT = ADCOUNT + 1
               DTMP = B1%TMP_G - B1%TMP_F
               IF (ABS(QNET) > 0._EB .AND. ABS(DTMP) <TWO_EPSILON_EB) DTMP=1._EB
               IF (PRESENT(WALL_INDEX)) THEN
                  B1%HEAT_TRANS_COEF = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC,SF,WALL_INDEX_IN=WALL_INDEX)
               ELSEIF (PRESENT(CFACE_INDEX)) THEN
                  B1%HEAT_TRANS_COEF = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC,SF,CFACE_INDEX_IN=CFACE_INDEX)
               ELSEIF (PRESENT(PARTICLE_INDEX)) THEN
                  B1%HEAT_TRANS_COEF = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC,SF,PARTICLE_INDEX_IN=PARTICLE_INDEX)
               ENDIF
               ! Use Ferrari's method
               IF (.NOT. RADIATION) THEN
                  IF (ABS(QNET) < TWO_EPSILON_EB) THEN
                     HTC_OLD = 0._EB
                     B1%TMP_F = B1%TMP_G
                  ELSE
                     HTC_OLD = B1%HEAT_TRANS_COEF
                     IF (ABS(B1%HEAT_TRANS_COEF) < TWO_EPSILON_EB) EXIT ADLOOP
                     B1%TMP_F = (-QNET + B1%HEAT_TRANS_COEF * B1%TMP_G)/B1%HEAT_TRANS_COEF
                  ENDIF
                  EXIT ADLOOP
               ENDIF
               BBB = B1%HEAT_TRANS_COEF / (B1%EMISSIVITY * SIGMA)
               CCC = -(B1%Q_RAD_IN+B1%HEAT_TRANS_COEF*B1%TMP_G-QNET)/(B1%EMISSIVITY * SIGMA)
               PPP = -CCC
               QQQ = -0.125_EB*BBB**2
               RRR = -0.5_EB*QQQ+SQRT(0.25_EB*QQQ**2+PPP**3/27._EB)
               UUU = RRR**ONTH
               IF (UUU < TWO_EPSILON_EB) THEN
                  YYY = -QQQ**ONTH
               ELSE
                  YYY = UUU-ONTH*PPP/UUU
               ENDIF
               WWW = SQRT(2._EB*YYY)
               B1%TMP_F = 0.5_EB*(-WWW+SQRT(-2._EB*(YYY-BBB/WWW)))
               HTC_OLD = 0.2_EB*HTC_OLD+0.8_EB*B1%HEAT_TRANS_COEF
               IF (ABS(HTC_OLD-B1%HEAT_TRANS_COEF) < 1.E-6_EB .OR. ADCOUNT > 20) EXIT ADLOOP
            ENDDO ADLOOP
            B1%HEAT_TRANS_COEF = HTC_OLD
            B1%Q_CON_F = B1%HEAT_TRANS_COEF*(B1%TMP_G-B1%TMP_F)
         ENDIF
      ENDIF
      IF (RADIATION) B1%Q_RAD_OUT = SIGMA*B1%EMISSIVITY*B1%TMP_F**4

   CASE (INTERPOLATED_BC) METHOD_OF_HEAT_TRANSFER  ! Only for EXTERNAL_WALL_CELLs

      EWC => EXTERNAL_WALL(WALL_INDEX)
      OM => OMESH(EWC%NOM)
      IF (PREDICTOR) THEN
         OM_RHOP => OM%RHOS
         OM_ZZP => OM%ZZS
      ELSE
         OM_RHOP => OM%RHO
         OM_ZZP => OM%ZZ
      ENDIF

      MM => MESHES(EWC%NOM)
      RHO_G_2     = B1%RHO_G ! first-order
      RHO_OTHER   = RHOP(BC%II,BC%JJ,BC%KK)
      RHO_OTHER_2 = RHOP(BC%II,BC%JJ,BC%KK)
      RHO_ZZ_OTHER(1:N_TOTAL_SCALARS) = ZZP(BC%II,BC%JJ,BC%KK,1:N_TOTAL_SCALARS)*RHO_OTHER

      ! Determine if there are 4 equally sized cells spanning the interpolated boundary

      SECOND_ORDER_INTERPOLATED_BOUNDARY = .FALSE.
      CC_SOLID_FLAG = .FALSE.
      IF (ABS(EWC%AREA_RATIO-1._EB)<0.01_EB) THEN
         IIO = EWC%IIO_MIN
         JJO = EWC%JJO_MIN
         KKO = EWC%KKO_MIN
         SELECT CASE(BC%IOR)
            CASE( 1) ; ICG = CELL_INDEX(BC%IIG+1,BC%JJG,BC%KKG) ; ICO = MM%CELL_INDEX(IIO-1,JJO,KKO)
            CASE(-1) ; ICG = CELL_INDEX(BC%IIG-1,BC%JJG,BC%KKG) ; ICO = MM%CELL_INDEX(IIO+1,JJO,KKO)
            CASE( 2) ; ICG = CELL_INDEX(BC%IIG,BC%JJG+1,BC%KKG) ; ICO = MM%CELL_INDEX(IIO,JJO-1,KKO)
            CASE(-2) ; ICG = CELL_INDEX(BC%IIG,BC%JJG-1,BC%KKG) ; ICO = MM%CELL_INDEX(IIO,JJO+1,KKO)
            CASE( 3) ; ICG = CELL_INDEX(BC%IIG,BC%JJG,BC%KKG+1) ; ICO = MM%CELL_INDEX(IIO,JJO,KKO-1)
            CASE(-3) ; ICG = CELL_INDEX(BC%IIG,BC%JJG,BC%KKG-1) ; ICO = MM%CELL_INDEX(IIO,JJO,KKO+1)
         END SELECT
         IF (CC_IBM) THEN ! Test if one of surrounding cells is CC_SOLID.
            IF(CCVAR(BC%IIG,BC%JJG,BC%KKG,CC_CGSC)==CC_SOLID .OR. CCVAR(BC%II,BC%JJ,BC%KK,CC_CGSC)==CC_SOLID) CC_SOLID_FLAG = .TRUE.
         ENDIF
         IF (.NOT.CELL(ICG)%SOLID.AND..NOT.MM%CELL(ICO)%SOLID.AND..NOT.CC_SOLID_FLAG) SECOND_ORDER_INTERPOLATED_BOUNDARY = .TRUE.
      ENDIF

      ! Density

      ATMOSPHERIC_INTERPOLATION = .FALSE.
      DDO = 1._EB
      KKO = EWC%KKO_MIN
      SELECT CASE(BC%IOR)
         CASE( 3) ; DDO = (DZ(BC%KK)+DZ(BC%KKG))/(MM%DZ(KKO)+DZ(BC%KKG))
         CASE(-3) ; DDO = (DZ(BC%KK)+DZ(BC%KKG))/(MM%DZ(KKO)+DZ(BC%KKG))
      END SELECT
      IF (USE_ATMOSPHERIC_INTERPOLATION .AND. STRATIFICATION .AND. ABS(DDO-1._EB)>0.01_EB .AND. ABS(BC%IOR)==3) &
         ATMOSPHERIC_INTERPOLATION = .TRUE.

      IF (ATMOSPHERIC_INTERPOLATION) THEN
         ! interp or extrap RHO_OTHER for jump in vertical grid resolution, linear in temperature to match heat flux in divg
         PBAR_G = PBAR_P(BC%KKG,B1%PRESSURE_ZONE)
         PBAR_OTHER = EVALUATE_RAMP(MM%ZC(EWC%KKO_MIN),I_RAMP_P0_Z)
         ZZ_GET(1:N_TRACKED_SPECIES) = MAX(0._EB,MIN(1._EB,ZZP(BC%II,BC%JJ,BC%KK,1:N_TOTAL_SCALARS)/RHOP(BC%II,BC%JJ,BC%KK)))
         CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET,RSUM(BC%II,BC%JJ,BC%KK))
         DENOM = PBAR_G/B1%RHO_G/RSUM_G + DDO*(PBAR_OTHER/RHO_OTHER/RSUM(BC%II,BC%JJ,BC%KK) - PBAR_G/B1%RHO_G/RSUM_G)
         RHOP(BC%II,BC%JJ,BC%KK) = PBAR_P(BC%KK,B1%PRESSURE_ZONE)/RSUM(BC%II,BC%JJ,BC%KK)/DENOM
      ENDIF

      SELECT CASE(BC%IOR)
         CASE( 1)
            IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN
               RHO_G_2 = RHOP(BC%IIG+1,BC%JJG,BC%KKG)
               RHO_OTHER_2 = OM_RHOP(IIO-1,JJO,KKO)
            ENDIF
            Z_TEMP(0:3,1,1) = (/RHO_OTHER_2,RHO_OTHER,B1%RHO_G,RHO_G_2/)
            U_TEMP(1,1,1) = UU(BC%II,BC%JJ,BC%KK)
            CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER)
            B1%RHO_F = F_TEMP(1,1,1)
         CASE(-1)
            IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN
               RHO_G_2 = RHOP(BC%IIG-1,BC%JJG,BC%KKG)
               RHO_OTHER_2 = OM_RHOP(IIO+1,JJO,KKO)
            ENDIF
            Z_TEMP(0:3,1,1) = (/RHO_G_2,B1%RHO_G,RHO_OTHER,RHO_OTHER_2/)
            U_TEMP(1,1,1) = UU(BC%II-1,BC%JJ,BC%KK)
            CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER)
            B1%RHO_F = F_TEMP(1,1,1)
         CASE( 2)
            IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN
               RHO_G_2 = RHOP(BC%IIG,BC%JJG+1,BC%KKG)
               RHO_OTHER_2 = OM_RHOP(IIO,JJO-1,KKO)
            ENDIF
            Z_TEMP(1,0:3,1) = (/RHO_OTHER_2,RHO_OTHER,B1%RHO_G,RHO_G_2/)
            U_TEMP(1,1,1) = VV(BC%II,BC%JJ,BC%KK)
            CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER)
            B1%RHO_F = F_TEMP(1,1,1)
         CASE(-2)
            IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN
               RHO_G_2 = RHOP(BC%IIG,BC%JJG-1,BC%KKG)
               RHO_OTHER_2 = OM_RHOP(IIO,JJO+1,KKO)
            ENDIF
            Z_TEMP(1,0:3,1) = (/RHO_G_2,B1%RHO_G,RHO_OTHER,RHO_OTHER_2/)
            U_TEMP(1,1,1) = VV(BC%II,BC%JJ-1,BC%KK)
            CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER)
            B1%RHO_F = F_TEMP(1,1,1)
         CASE( 3)
            IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN
               RHO_G_2 = RHOP(BC%IIG,BC%JJG,BC%KKG+1)
               RHO_OTHER_2 = OM_RHOP(IIO,JJO,KKO-1)
            ENDIF
            Z_TEMP(1,1,0:3) = (/RHO_OTHER_2,RHO_OTHER,B1%RHO_G,RHO_G_2/)
            U_TEMP(1,1,1) = WW(BC%II,BC%JJ,BC%KK)
            CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER)
            B1%RHO_F = F_TEMP(1,1,1)
         CASE(-3)
            IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN
               RHO_G_2 = RHOP(BC%IIG,BC%JJG,BC%KKG-1)
               RHO_OTHER_2 = OM_RHOP(IIO,JJO,KKO+1)
            ENDIF
            Z_TEMP(1,1,0:3) = (/RHO_G_2,B1%RHO_G,RHO_OTHER,RHO_OTHER_2/)
            U_TEMP(1,1,1) = WW(BC%II,BC%JJ,BC%KK-1)
            CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER)
            B1%RHO_F = F_TEMP(1,1,1)
      END SELECT

      ! Species and temperature

      SINGLE_SPEC_IF: IF (N_TOTAL_SCALARS > 1) THEN
         SPECIES_LOOP: DO N=1,N_TOTAL_SCALARS

            RHO_ZZ_G = B1%RHO_G*B1%ZZ_G(N)
            RHO_ZZ_G_2 = RHO_ZZ_G ! first-order (default)
            RHO_ZZ_OTHER_2 = RHO_ZZ_OTHER(N)

            SELECT CASE(BC%IOR)
               CASE( 1)
                  IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN
                     RHO_ZZ_G_2 = RHOP(BC%IIG+1,BC%JJG,BC%KKG)*ZZP(BC%IIG+1,BC%JJG,BC%KKG,N)
                     RHO_ZZ_OTHER_2 = OM_RHOP(IIO-1,JJO,KKO)*OM_ZZP(IIO-1,JJO,KKO,N)
                  ENDIF
                  Z_TEMP(0:3,1,1) = (/RHO_ZZ_OTHER_2,RHO_ZZ_OTHER(N),RHO_ZZ_G,RHO_ZZ_G_2/)
                  U_TEMP(1,1,1) = UU(BC%II,BC%JJ,BC%KK)
                  CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER)
                  RHO_ZZ_F(N) = F_TEMP(1,1,1)
                  PBAR_F = PBAR_P(BC%KKG,B1%PRESSURE_ZONE)
               CASE(-1)
                  IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN
                     RHO_ZZ_G_2 = RHOP(BC%IIG-1,BC%JJG,BC%KKG)*ZZP(BC%IIG-1,BC%JJG,BC%KKG,N)
                     RHO_ZZ_OTHER_2 = OM_RHOP(IIO+1,JJO,KKO)*OM_ZZP(IIO+1,JJO,KKO,N)
                  ENDIF
                  Z_TEMP(0:3,1,1) = (/RHO_ZZ_G_2,RHO_ZZ_G,RHO_ZZ_OTHER(N),RHO_ZZ_OTHER_2/)
                  U_TEMP(1,1,1) = UU(BC%II-1,BC%JJ,BC%KK)
                  CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,1,I_FLUX_LIMITER)
                  RHO_ZZ_F(N) = F_TEMP(1,1,1)
                  PBAR_F = PBAR_P(BC%KKG,B1%PRESSURE_ZONE)
               CASE( 2)
                  IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN
                     RHO_ZZ_G_2 = RHOP(BC%IIG,BC%JJG+1,BC%KKG)*ZZP(BC%IIG,BC%JJG+1,BC%KKG,N)
                     RHO_ZZ_OTHER_2 = OM_RHOP(IIO,JJO-1,KKO)*OM_ZZP(IIO,JJO-1,KKO,N)
                  ENDIF
                  Z_TEMP(1,0:3,1) = (/RHO_ZZ_OTHER_2,RHO_ZZ_OTHER(N),RHO_ZZ_G,RHO_ZZ_G_2/)
                  U_TEMP(1,1,1) = VV(BC%II,BC%JJ,BC%KK)
                  CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER)
                  RHO_ZZ_F(N) = F_TEMP(1,1,1)
                  PBAR_F = PBAR_P(BC%KKG,B1%PRESSURE_ZONE)
               CASE(-2)
                  IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN
                     RHO_ZZ_G_2 = RHOP(BC%IIG,BC%JJG-1,BC%KKG)*ZZP(BC%IIG,BC%JJG-1,BC%KKG,N)
                     RHO_ZZ_OTHER_2 = OM_RHOP(IIO,JJO+1,KKO)*OM_ZZP(IIO,JJO+1,KKO,N)
                  ENDIF
                  Z_TEMP(1,0:3,1) = (/RHO_ZZ_G_2,RHO_ZZ_G,RHO_ZZ_OTHER(N),RHO_ZZ_OTHER_2/)
                  U_TEMP(1,1,1) = VV(BC%II,BC%JJ-1,BC%KK)
                  CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,2,I_FLUX_LIMITER)
                  RHO_ZZ_F(N) = F_TEMP(1,1,1)
                  PBAR_F = PBAR_P(BC%KKG,B1%PRESSURE_ZONE)
               CASE( 3)
                  IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN
                     RHO_ZZ_G_2 = RHOP(BC%IIG,BC%JJG,BC%KKG+1)*ZZP(BC%IIG,BC%JJG,BC%KKG+1,N)
                     RHO_ZZ_OTHER_2 = OM_RHOP(IIO,JJO,KKO-1)*OM_ZZP(IIO,JJO,KKO-1,N)
                  ENDIF
                  Z_TEMP(1,1,0:3) = (/RHO_ZZ_OTHER_2,RHO_ZZ_OTHER(N),RHO_ZZ_G,RHO_ZZ_G_2/)
                  U_TEMP(1,1,1) = WW(BC%II,BC%JJ,BC%KK)
                  CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER)
                  RHO_ZZ_F(N) = F_TEMP(1,1,1)
                  PBAR_F = (PBAR_P(BC%KK,B1%PRESSURE_ZONE)*DZ(BC%KKG) + PBAR_P(BC%KKG,B1%PRESSURE_ZONE)*DZ(BC%KK)) / &
                           (DZ(BC%KK)+DZ(BC%KKG))
               CASE(-3)
                  IF (SECOND_ORDER_INTERPOLATED_BOUNDARY) THEN
                     RHO_ZZ_G_2 = RHOP(BC%IIG,BC%JJG,BC%KKG-1)*ZZP(BC%IIG,BC%JJG,BC%KKG-1,N)
                     RHO_ZZ_OTHER_2 = OM_RHOP(IIO,JJO,KKO+1)*OM_ZZP(IIO,JJO,KKO+1,N)
                  ENDIF
                  Z_TEMP(1,1,0:3) = (/RHO_ZZ_G_2,RHO_ZZ_G,RHO_ZZ_OTHER(N),RHO_ZZ_OTHER_2/)
                  U_TEMP(1,1,1) = WW(BC%II,BC%JJ,BC%KK-1)
                  CALL GET_SCALAR_FACE_VALUE(U_TEMP,Z_TEMP,F_TEMP,1,1,1,1,1,1,3,I_FLUX_LIMITER)
                  RHO_ZZ_F(N) = F_TEMP(1,1,1)
                  PBAR_F = (PBAR_P(BC%KK,B1%PRESSURE_ZONE)*DZ(BC%KKG) + PBAR_P(BC%KKG,B1%PRESSURE_ZONE)*DZ(BC%KK)) / &
                           (DZ(BC%KK)+DZ(BC%KKG))
            END SELECT
         ENDDO SPECIES_LOOP

         ! face value of temperature

         IF (ATMOSPHERIC_INTERPOLATION) THEN
            B1%TMP_F = (TMP(BC%II,BC%JJ,BC%KK)*DZ(BC%KKG) + TMP(BC%IIG,BC%JJG,BC%KKG)*DZ(BC%KK)) / (DZ(BC%KK)+DZ(BC%KKG))
            B1%ZZ_F(1:N_TOTAL_SCALARS) = (ZZP(BC%II,BC%JJ,BC%KK,1:N_TOTAL_SCALARS)*DZ(BC%KKG) + &
                                          ZZP(BC%IIG,BC%JJG,BC%KKG,1:N_TOTAL_SCALARS)*DZ(BC%KK)) / (DZ(BC%KK)+DZ(BC%KKG))
            ZZ_GET(1:N_TRACKED_SPECIES) = B1%ZZ_F(1:N_TRACKED_SPECIES)
            CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET,RSUM_F)
            B1%RHO_F = PBAR_F/(RSUM_F*B1%TMP_F)
         ELSE
            B1%ZZ_F(1:N_TOTAL_SCALARS) = MAX(0._EB,MIN(1._EB,RHO_ZZ_F(1:N_TOTAL_SCALARS)/B1%RHO_F))
            ZZ_GET(1:N_TRACKED_SPECIES) = B1%ZZ_F(1:N_TRACKED_SPECIES)
            CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET,RSUM_F)
            B1%TMP_F = PBAR_F/(RSUM_F*B1%RHO_F)
         ENDIF

         ! flux match species diffusive flux at interpolated boundaries with mesh refinement
         COARSE_MESH_IF: IF (EWC%NIC>1) THEN
            ! we are on coarse mesh gas cell (G) and need to average fluxes from the fine mesh (OTHER)
            OM_MUP => OM%MU
            SPECIES_LOOP_2: DO N=1,N_TOTAL_SCALARS
               SELECT CASE(SIM_MODE)
                  CASE(DNS_MODE,LES_MODE)
                     D_Z_N = D_Z(:,N)
                     CALL INTERPOLATE1D_UNIFORM(LBOUND(D_Z_N,1),D_Z_N,B1%TMP_G,D_Z_G)
                     IF (SIM_MODE==LES_MODE) CALL GET_VISCOSITY(ZZ_GET,MU_DNS_G,B1%TMP_G)
               END SELECT
               RHO_D_DZDN_OTHER = 0._EB
               KKO_LOOP: DO KKO=EWC%KKO_MIN,EWC%KKO_MAX
                  JJO_LOOP: DO JJO=EWC%JJO_MIN,EWC%JJO_MAX
                     IIO_LOOP: DO IIO=EWC%IIO_MIN,EWC%IIO_MAX
                        MU_OTHER = OM_MUP(IIO,JJO,KKO)
                        RHO_OTHER = OM_RHOP(IIO,JJO,KKO)
                        MODE_SELECT: SELECT CASE(SIM_MODE)
                           CASE DEFAULT
                              RHO_D = MAX(0._EB, 0.5_EB*(MU_OTHER+MU_G) )*RSC
                           CASE(DNS_MODE,LES_MODE)
                              D_Z_N = D_Z(:,N)
                              ZZ_GET(1:N_TRACKED_SPECIES) = OM_ZZP(IIO,JJO,KKO,1:N_TRACKED_SPECIES)
                              CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET,RSUM_OTHER)
                              PBAR_OTHER = EVALUATE_RAMP(MM%ZC(KKO),I_RAMP_P0_Z)
                              TMP_OTHER = PBAR_OTHER/(RSUM_OTHER*RHO_OTHER)
                              CALL INTERPOLATE1D_UNIFORM(LBOUND(D_Z_N,1),D_Z_N,TMP_OTHER,D_Z_OTHER)
                              RHO_D_TURB = 0._EB
                              IF (SIM_MODE==LES_MODE) THEN
                                 CALL GET_VISCOSITY(ZZ_GET,MU_DNS_OTHER,TMP_OTHER)
                                 RSC_LOC = RSC
                                 IF (SPECIES_MIXTURE(N)%SC_T_USER>TWO_EPSILON_EB) RSC_LOC=1._EB/SPECIES_MIXTURE(N)%SC_T_USER
                                 RHO_D_TURB = 0.5_EB*(MU_OTHER-MU_DNS_OTHER + MU_G-MU_DNS_G)*RSC_LOC
                              ENDIF
                              RHO_D = 0.5_EB*( RHO_OTHER*D_Z_OTHER + B1%RHO_G*D_Z_G ) + RHO_D_TURB
                        END SELECT MODE_SELECT
                        SELECT CASE(BC%IOR)
                           CASE( 1)
                              ARO = MM%DY(JJO)*MM%DZ(KKO)/(DY(BC%JJ)*DZ(BC%KK))
                              RHO_D_DZDN = RHO_D*(ZZP(BC%IIG,BC%JJG,BC%KKG,N)-OM_ZZP(IIO,JJO,KKO,N))*MM%RDXN(IIO)
                           CASE(-1)
                              ARO = MM%DY(JJO)*MM%DZ(KKO)/(DY(BC%JJ)*DZ(BC%KK))
                              RHO_D_DZDN = RHO_D*(ZZP(BC%IIG,BC%JJG,BC%KKG,N)-OM_ZZP(IIO,JJO,KKO,N))*MM%RDXN(IIO-1)
                           CASE( 2)
                              ARO = MM%DX(IIO)*MM%DZ(KKO)/(DX(BC%II)*DZ(BC%KK))
                              RHO_D_DZDN = RHO_D*(ZZP(BC%IIG,BC%JJG,BC%KKG,N)-OM_ZZP(IIO,JJO,KKO,N))*MM%RDYN(JJO)
                           CASE(-2)
                              ARO = MM%DX(IIO)*MM%DZ(KKO)/(DX(BC%II)*DZ(BC%KK))
                              RHO_D_DZDN = RHO_D*(ZZP(BC%IIG,BC%JJG,BC%KKG,N)-OM_ZZP(IIO,JJO,KKO,N))*MM%RDYN(JJO-1)
                           CASE( 3)
                              ARO = MM%DX(IIO)*MM%DY(JJO)/(DX(BC%II)*DY(BC%JJ))
                              RHO_D_DZDN = RHO_D*(ZZP(BC%IIG,BC%JJG,BC%KKG,N)-OM_ZZP(IIO,JJO,KKO,N))*MM%RDZN(KKO)
                           CASE(-3)
                              ARO = MM%DX(IIO)*MM%DY(JJO)/(DX(BC%II)*DY(BC%JJ))
                              RHO_D_DZDN = RHO_D*(ZZP(BC%IIG,BC%JJG,BC%KKG,N)-OM_ZZP(IIO,JJO,KKO,N))*MM%RDZN(KKO-1)
                        END SELECT
                        ! average multiple face values
                        RHO_D_DZDN_OTHER = RHO_D_DZDN_OTHER + ARO*RHO_D_DZDN
                     ENDDO IIO_LOOP
                  ENDDO JJO_LOOP
               ENDDO KKO_LOOP
               ! store for use in divg
               B1%RHO_D_DZDN_F(N) =  RHO_D_DZDN_OTHER
            ENDDO SPECIES_LOOP_2
            IF (SIM_MODE==DNS_MODE .OR. SIM_MODE==LES_MODE) THEN
               N=MAXLOC(B1%ZZ_F(1:N_TRACKED_SPECIES),1)
               B1%RHO_D_DZDN_F(N) = -(SUM(B1%RHO_D_DZDN_F(1:N_TRACKED_SPECIES))-B1%RHO_D_DZDN_F(N))
            ENDIF
         ENDIF COARSE_MESH_IF

      ELSE SINGLE_SPEC_IF

         B1%ZZ_F(1) = 1._EB
         TMP(BC%II,BC%JJ,BC%KK) = PBAR_P(BC%KK,B1%PRESSURE_ZONE)/(RSUM0*RHOP(BC%II,BC%JJ,BC%KK))
         SELECT CASE(BC%IOR)
            CASE DEFAULT
               PBAR_F = PBAR_P(BC%KKG,B1%PRESSURE_ZONE)
            CASE (-3,3)
               PBAR_F = (PBAR_P(BC%KK,B1%PRESSURE_ZONE)*DZ(BC%KKG) + PBAR_P(BC%KKG,B1%PRESSURE_ZONE)*DZ(BC%KK)) / &
                        (DZ(BC%KK)+DZ(BC%KKG))
         END SELECT
         IF (ATMOSPHERIC_INTERPOLATION) THEN
            B1%TMP_F = (TMP(BC%II,BC%JJ,BC%KK)*DZ(BC%KKG) + TMP(BC%IIG,BC%JJG,BC%KKG)*DZ(BC%KK)) / (DZ(BC%KK)+DZ(BC%KKG))
            B1%RHO_F = PBAR_F/(RSUM0*B1%TMP_F)
         ELSE
            B1%TMP_F = PBAR_F/(RSUM0*B1%RHO_F)
         ENDIF
      ENDIF SINGLE_SPEC_IF

      B1%Q_CON_F = 0._EB ! no convective heat transfer at interpolated boundary
      IF (RADIATION) B1%Q_RAD_OUT = SIGMA*B1%EMISSIVITY*B1%TMP_F**4

END SELECT METHOD_OF_HEAT_TRANSFER

END SUBROUTINE SURFACE_HEAT_TRANSFER


!> \brief Calculate the diffusion coefficient at the boundary
!> \param B1 Pointer to boundary properties
!> \param BC Pointer to boundary coordinates
!> \param WALL_INDEX Optional WALL cell index
!> \param CFACE_INDEX Optional immersed boundary (CFACE) index
!> \details Calculate the diffusion coefficient, rho*D (kg/m/s)

SUBROUTINE CALCULATE_RHO_D_F(B1,BC,WALL_INDEX,CFACE_INDEX)

INTEGER, INTENT(IN), OPTIONAL :: WALL_INDEX,CFACE_INDEX
REAL(EB) :: MU_G
TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1
TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC
INTEGER :: N,ITMP
REAL(EB) :: RSC_LOC

IF (PRESENT(WALL_INDEX)) THEN
   B1%RHO_G = RHO(BC%IIG,BC%JJG,BC%KKG)
   MU_G = MU(BC%IIG,BC%JJG,BC%KKG)
ELSEIF (PRESENT(CFACE_INDEX)) THEN
   MU_G  = CFACE(CFACE_INDEX)%MU_G
ENDIF

SELECT CASE(SIM_MODE)
   CASE DEFAULT
      DO N=1,N_TRACKED_SPECIES
         B1%RHO_D_F(N) = MU_G*RSC*B1%RHO_F/B1%RHO_G
      ENDDO
   CASE (LES_MODE)
      ITMP = MIN(I_MAX_TEMP-1,NINT(B1%TMP_F))
      DO N=1,N_TRACKED_SPECIES
         RSC_LOC = RSC
         IF (SPECIES_MIXTURE(N)%SC_T_USER>TWO_EPSILON_EB) RSC_LOC=1._EB/SPECIES_MIXTURE(N)%SC_T_USER
         B1%RHO_D_F(N) = B1%RHO_F*( D_Z(ITMP,N) + (MU_G-MU_DNS(BC%IIG,BC%JJG,BC%KKG))/B1%RHO_G*RSC_LOC )
      ENDDO
   CASE (DNS_MODE)
      ITMP = MIN(I_MAX_TEMP-1,NINT(B1%TMP_F))
      DO N=1,N_TRACKED_SPECIES
         B1%RHO_D_F(N) = B1%RHO_F*D_Z(ITMP,N)
      ENDDO
END SELECT

END SUBROUTINE CALCULATE_RHO_D_F


!> \brief Calculate the species mass fractions, ZZ, at the boundary
!> \param T Current time (s)
!> \param DT Current time step (s)
!> \param WALL_INDEX Optional WALL cell index
!> \param CFACE_INDEX Optional immersed boundary (CFACE) index
!> \param PARTICLE_INDEX Optional particle index
!> \details Calculate the diffusion coefficient, rho*D (kg/m/s)

SUBROUTINE CALCULATE_ZZ_F(T,DT,WALL_INDEX,CFACE_INDEX,PARTICLE_INDEX)

USE HVAC_ROUTINES, ONLY : DUCT_MF
USE PHYSICAL_FUNCTIONS, ONLY: GET_SPECIFIC_HEAT,GET_SPECIFIC_GAS_CONSTANT, GET_REALIZABLE_MF, Q_REF_FIT
USE MATH_FUNCTIONS, ONLY : EVALUATE_RAMP, BOX_MULLER, INTERPOLATE1D_UNIFORM
REAL(EB), INTENT(IN) :: T,DT
REAL(EB) :: UN,DD,MFT,TSI,RSUM_F,MPUA_SUM,RHO_F_PREVIOUS,RN1,RN2,MFT_UNIFORM,Q_NEW(MAX_QDOTPP_REF)
REAL(EB) :: T_SCALE(MAX_QDOTPP_REF),QDOTPP_REF(MAX_QDOTPP_REF),QDOTPP_T(MAX_QDOTPP_REF), &
            QDOTPP,QDOTPP1,QDOTPP2,DT_SPYRO(MAX_QDOTPP_REF),CP,H_G,MW_RATIO
REAL(EB) :: RVC,M_DOT_PPP_SINGLE,ZZ_GET(1:N_TRACKED_SPECIES),DENOM
INTEGER :: N,NS,IDX1,IDX2,NQ,ITER,IIO,JJO,KKO,OBST_INDEX,OTHER_MESH_OBST_INDEX,LL,SPECIES_BC_INDEX,IC,ICG
INTEGER, INTENT(IN), OPTIONAL :: WALL_INDEX,CFACE_INDEX,PARTICLE_INDEX
TYPE(RAMPS_TYPE), POINTER :: RP,RP_E2T,RP_REF
TYPE(EXTERNAL_WALL_TYPE), POINTER :: EWC
TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1
TYPE(BOUNDARY_PROP2_TYPE), POINTER :: B2
TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC
TYPE(SURFACE_TYPE), POINTER :: SF
TYPE(LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC
TYPE(LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP
TYPE(CFACE_TYPE), POINTER :: CFA
TYPE(WALL_TYPE), POINTER :: WC

IF (PRESENT(WALL_INDEX)) THEN
   WC => WALL(WALL_INDEX)
   B1 => BOUNDARY_PROP1(WC%B1_INDEX)
   B2 => BOUNDARY_PROP2(WC%B2_INDEX)
   BC => BOUNDARY_COORD(WC%BC_INDEX)
   SF => SURFACE(WC%SURF_INDEX)
   IC  = CELL_INDEX(BC%II,BC%JJ,BC%KK)
   ICG = CELL_INDEX(BC%IIG,BC%JJG,BC%KKG)
   OBST_INDEX = WALL(WALL_INDEX)%OBST_INDEX
ELSEIF (PRESENT(CFACE_INDEX)) THEN
   CFA => CFACE(CFACE_INDEX)
   B1 => BOUNDARY_PROP1(CFA%B1_INDEX)
   B2 => BOUNDARY_PROP2(CFA%B2_INDEX)
   BC => BOUNDARY_COORD(CFA%BC_INDEX)
   SF => SURFACE(CFA%SURF_INDEX)
   IC  = 0
   ICG = 0
   OBST_INDEX = 0
ELSEIF (PRESENT(PARTICLE_INDEX)) THEN
   LP  => LAGRANGIAN_PARTICLE(PARTICLE_INDEX)
   LPC => LAGRANGIAN_PARTICLE_CLASS(LP%CLASS_INDEX)
   SF => SURFACE(LPC%SURF_INDEX)
   B1 => BOUNDARY_PROP1(LP%B1_INDEX)
   BC => BOUNDARY_COORD(LP%BC_INDEX)
   IC  = CELL_INDEX(BC%IIG,BC%JJG,BC%KKG)
   OBST_INDEX = 0
ENDIF

! Special cases for N_TRACKED_SPECIES==1

IF (N_TRACKED_SPECIES==1) THEN

   IF ( B1%NODE_INDEX < 0 .AND. .NOT.SF%SPECIES_BC_INDEX==SPECIFIED_MASS_FLUX ) THEN
      B1%ZZ_F(1) = 1._EB
      RETURN
   ENDIF

   IF ( SF%SPECIES_BC_INDEX==SPECIFIED_MASS_FLUX .AND. ABS(SF%MASS_FLUX(1))<=TWO_EPSILON_EB ) THEN
      B1%ZZ_F(1) = 1._EB
      RETURN
   ENDIF

ENDIF

! Check if suppression by water is to be applied and sum water on surface

IF (CORRECTOR .AND. SF%E_COEFFICIENT>0._EB .AND. I_WATER>0 .AND. (PRESENT(WALL_INDEX).OR.PRESENT(CFACE_INDEX))) THEN
   IF (SPECIES_MIXTURE(I_WATER)%EVAPORATING) THEN
      MPUA_SUM = 0._EB
      DO N=1,N_LAGRANGIAN_CLASSES
         LPC=>LAGRANGIAN_PARTICLE_CLASS(N)
         IF (LPC%Z_INDEX==I_WATER) MPUA_SUM = MPUA_SUM + B2%LP_MPUA(LPC%ARRAY_INDEX)
      ENDDO
      B2%K_SUPPRESSION = B2%K_SUPPRESSION + SF%E_COEFFICIENT*MPUA_SUM*DT
   ENDIF
ENDIF

! Adjust SPECIES_BC_INDEX for HVAC

SPECIES_BC_INDEX = SF%SPECIES_BC_INDEX

IF (B1%NODE_INDEX > 0) THEN
   IF (-DUCTNODE(B1%NODE_INDEX)%DIR(1)*DUCT_MF(DUCTNODE(B1%NODE_INDEX)%DUCT_INDEX(1))>=0._EB) THEN
      SPECIES_BC_INDEX = SPECIFIED_MASS_FRACTION
   ELSE
      SPECIES_BC_INDEX = SPECIFIED_MASS_FLUX
   ENDIF
ENDIF

! Apply the different species boundary conditions to non-thermally thick solids

METHOD_OF_MASS_TRANSFER: SELECT CASE(SPECIES_BC_INDEX)

   CASE (INFLOW_OUTFLOW_MASS_FLUX) METHOD_OF_MASS_TRANSFER

      ! OPEN boundary species BC is done in THERMAL_BC under INFLOW_OUTFLOW

   CASE (NO_MASS_FLUX) METHOD_OF_MASS_TRANSFER

      B1%ZZ_F(1:N_TRACKED_SPECIES) = B1%ZZ_G(1:N_TRACKED_SPECIES)

   CASE (SPECIFIED_MASS_FRACTION) METHOD_OF_MASS_TRANSFER

      IF (ABS(B1%T_IGN-T_BEGIN)<SPACING(B1%T_IGN) .AND. ANY(SF%RAMP(1:N_TRACKED_SPECIES)%INDEX>=1)) THEN
         IF (PREDICTOR) TSI = T + DT
         IF (CORRECTOR) TSI = T
      ELSE
         IF (PREDICTOR) TSI = T + DT - B1%T_IGN
         IF (CORRECTOR) TSI = T      - B1%T_IGN
      ENDIF

      IF (B1%U_NORMAL_S<0._EB) THEN  ! If there is a non-zero velocity into the domain, assign appropriate species
                                        ! mass fractions to the face
         DO N=2,N_TRACKED_SPECIES
            ZZ_GET(N) = SPECIES_MIXTURE(N)%ZZ0 + EVALUATE_RAMP(TSI,SF%RAMP(N)%INDEX,TAU=SF%RAMP(N)%TAU)* &
                           (SF%MASS_FRACTION(N)-SPECIES_MIXTURE(N)%ZZ0)
         ENDDO
         ZZ_GET(1) = 1._EB-SUM(ZZ_GET(2:N_TRACKED_SPECIES))
         CALL GET_REALIZABLE_MF(ZZ_GET)
         B1%ZZ_F = ZZ_GET
      ELSE
         B1%ZZ_F(1:N_TRACKED_SPECIES) = B1%ZZ_G(1:N_TRACKED_SPECIES)
      ENDIF

      IF (PERIODIC_TEST==12 .AND. (TRIM(SF%ID)=='inlet')) THEN
         B1%ZZ_F(2) = 1._EB
         B1%ZZ_F(1) = 0._EB
      ENDIF
      IF (PERIODIC_TEST==13 .AND. (TRIM(SF%ID)=='inlet')) THEN
         B1%ZZ_F(2) = 0.5_EB*(1._EB + COS(4._EB*PI*XC(BC%II)))
         B1%ZZ_F(1) = 1._EB - B1%ZZ_F(2)
      ENDIF

      ! reconstruct species mass flux at the surface for output (some terms are lagged)

      SELECT CASE(BC%IOR)
         CASE( 1); UN = UU(BC%IIG-1,BC%JJG  ,BC%KKG  )
         CASE(-1); UN = UU(BC%IIG  ,BC%JJG  ,BC%KKG  )
         CASE( 2); UN = VV(BC%IIG  ,BC%JJG-1,BC%KKG  )
         CASE(-2); UN = VV(BC%IIG  ,BC%JJG  ,BC%KKG  )
         CASE( 3); UN = WW(BC%IIG  ,BC%JJG  ,BC%KKG-1)
         CASE(-3); UN = WW(BC%IIG  ,BC%JJG  ,BC%KKG  )
      END SELECT
      DO N=1,N_TRACKED_SPECIES
         B1%M_DOT_G_PP_ADJUST(N) = SIGN(1._EB,REAL(BC%IOR,EB))*( B1%RHO_F*B1%ZZ_F(N)*UN - B1%RHO_D_DZDN_F(N) )
         B1%M_DOT_G_PP_ACTUAL(N) = B1%M_DOT_G_PP_ADJUST(N)
      ENDDO

   CASE (SPECIFIED_MASS_FLUX) METHOD_OF_MASS_TRANSFER

      ! (S-Pyro) Calculate smoothed incident heat flux if cone scaling is applied

      IF (SF%N_QDOTPP_REF > 0 .AND. B1%T_IGN <=T .AND. PREDICTOR) THEN
         DT_SPYRO(1:SF%N_THICK_REF) = DT * SF%SPYRO_TH_FACTOR(1:SF%N_THICK_REF)
         TSI = MIN(T-B1%T_IGN+DT, SF%REFERENCE_HEAT_FLUX_TIME_INTERVAL+DT)
         IF (SOLID_PHASE_ONLY .AND. .NOT. SF%INERT_Q_REF) THEN
            B1%Q_IN_SMOOTH = (B1%Q_IN_SMOOTH *(TSI-DT) + &
               DT*Q_REF_FIT(SUM(B1%M_DOT_G_PP_ACTUAL)*SF%HOC_EFF,SF%HOC_EFF,SF%Y_S_EFF,B1%Q_RAD_IN/B1%EMISSIVITY))/TSI
         ELSE
            IF (B1%EMISSIVITY > 0._EB) THEN
               B1%Q_IN_SMOOTH = (B1%Q_IN_SMOOTH *(TSI-DT) + DT*(B1%Q_CON_F+B1%Q_RAD_IN/B1%EMISSIVITY))/TSI
            ELSE
               B1%Q_IN_SMOOTH = (B1%Q_IN_SMOOTH *(TSI-DT) + DT*(B1%Q_CON_F+B1%Q_RAD_IN))/TSI
            ENDIF
         ENDIF
         B1%Q_IN_SMOOTH = MIN(MAX(SF%MINIMUM_SCALING_HEAT_FLUX,B1%Q_IN_SMOOTH),SF%MAXIMUM_SCALING_HEAT_FLUX)
         B1%Q_IN_SMOOTH_INT(1:SF%N_THICK_REF) = B1%Q_IN_SMOOTH_INT(1:SF%N_THICK_REF) + DT_SPYRO(1:SF%N_THICK_REF) * B1%Q_IN_SMOOTH
      ENDIF

      ! If the current time is before the "activation" time, T_IGN, apply simple BCs and get out

      IF (T < B1%T_IGN .OR. B1%T_IGN+B1%BURN_DURATION<T .OR. INITIALIZATION_PHASE) THEN
         B1%ZZ_F(1:N_TRACKED_SPECIES) = B1%ZZ_G(1:N_TRACKED_SPECIES)
         IF (PREDICTOR) B1%U_NORMAL_S = 0._EB
         IF (CORRECTOR) B1%U_NORMAL  = 0._EB
         B1%M_DOT_G_PP_ADJUST(1:N_TRACKED_SPECIES) = 0._EB
         B1%M_DOT_G_PP_ACTUAL(1:N_TRACKED_SPECIES) = 0._EB
         B1%M_DOT_PART_ACTUAL = 0._EB
         RETURN
      ENDIF

      ! Zero out the running counter of Mass Flux Total (MFT)

      MFT = 0._EB

      ! (S-Pyro) If the user has specified the burning rate, evaluate the ramp and other related parameters

      IF (SF%N_QDOTPP_REF > 0 .AND. N_REACTIONS>=1 .AND. PREDICTOR .AND. B1%T_IGN <=T) THEN
         DO N = 1,SF%N_THICK_REF
            DO IDX1 = 1, SF%THICK2QREF(N,0)
               NQ = SF%THICK2QREF(N,IDX1)
               RP     => RAMPS(SF%HRRPUA_INT_INDEX(NQ))
               RP_REF => RAMPS(SF%QREF_INDEX(NQ))
               RP_E2T => RAMPS(SF%E2T_INDEX(NQ))
               CALL INTERPOLATE1D_UNIFORM(0,RP_E2T%INTERPOLATED_DATA(:),B1%Q_IN_SMOOTH_INT(N)*RP_E2T%RDT,T_SCALE(NQ))
               CALL INTERPOLATE1D_UNIFORM(0,RP_REF%INTERPOLATED_DATA(:),T_SCALE(NQ)*RP_REF%RDT,QDOTPP_REF(NQ))
               IF (B1%Q_IN_SMOOTH_INT(N) >= RP_E2T%T_MAX) THEN
                  RP => RAMPS(SF%RAMP(TIME_HEAT-NQ+1)%INDEX)
                  Q_NEW(NQ) = B1%QDOTPP_INT(NQ) + &
                                DT_SPYRO(N)*RP%INTERPOLATED_DATA(RP%NUMBER_INTERPOLATION_POINTS)*B1%Q_IN_SMOOTH/QDOTPP_REF(NQ)
               ELSE
                  CALL INTERPOLATE1D_UNIFORM(0,RP%INTERPOLATED_DATA(:),T_SCALE(NQ)*RP%RDT,Q_NEW(NQ))
               ENDIF
            ENDDO
         ENDDO
         DO N = 1,SF%N_THICK_REF
            IDX1 = SF%THICK2QREF(N,1)
            IDX2 = SF%THICK2QREF(N,SF%THICK2QREF(N,0))
            IF (B1%Q_IN_SMOOTH <= QDOTPP_REF(IDX1)) THEN
               QDOTPP_T(N) = MAX(0._EB,(Q_NEW(IDX1) - B1%QDOTPP_INT(IDX1))/DT_SPYRO(N))
            ELSEIF (B1%Q_IN_SMOOTH > QDOTPP_REF(IDX2)) THEN
               QDOTPP_T(N) = MAX(0._EB,(Q_NEW(IDX2) - B1%QDOTPP_INT(IDX2))/DT_SPYRO(N))
            ELSE
               DO NQ=IDX1+1,IDX2
                  IF (B1%Q_IN_SMOOTH <= QDOTPP_REF(NQ)) THEN
                     QDOTPP1 = MAX(0._EB,(Q_NEW(NQ-1)- B1%QDOTPP_INT(NQ-1))/DT_SPYRO(N))
                     QDOTPP2 = MAX(0._EB,(Q_NEW(NQ)  - B1%QDOTPP_INT(NQ))/DT_SPYRO(N))
                     QDOTPP_T(N) = QDOTPP1 + (QDOTPP2 - QDOTPP1)*(B1%Q_IN_SMOOTH - QDOTPP_REF(NQ-1))/&
                                   (QDOTPP_REF(NQ)-QDOTPP_REF(NQ-1))
                     EXIT
                  ENDIF
               ENDDO
            ENDIF
         ENDDO
         IF (SF%SPYRO_TH_FACTOR(1) >= 1._EB) THEN
            QDOTPP = QDOTPP_T(1)
         ELSE
            IF (SF%SPYRO_TH_FACTOR(SF%N_THICK_REF) <= 1._EB) THEN
               QDOTPP = QDOTPP_T(SF%N_THICK_REF)
            ELSE
               DO N=2,SF%N_THICK_REF
                  IF (SF%SPYRO_TH_FACTOR(N) > 1._EB) THEN
                     QDOTPP = QDOTPP_T(N-1) + (QDOTPP_T(N) - QDOTPP_T(N-1)) * (1._EB-SF%SPYRO_TH_FACTOR(N-1)) / &
                                              (SF%SPYRO_TH_FACTOR(N) - SF%SPYRO_TH_FACTOR(N-1))
                     EXIT
                  ENDIF
               ENDDO
            ENDIF
         ENDIF
         B1%QDOTPP_INT(1:SF%N_QDOTPP_REF) = Q_NEW(1:SF%N_QDOTPP_REF)
      ENDIF

      ! Determine the mass flux total (MFT) at solid obstruction

      SUM_MASSFLUX_LOOP: DO N=1,N_TRACKED_SPECIES
         IF (ABS(SF%MASS_FLUX(N)) > TWO_EPSILON_EB) THEN  ! Use user-specified ramp-up of mass flux
            IF (ABS(B1%T_IGN-T_BEGIN) < SPACING(B1%T_IGN) .AND. SF%RAMP(N)%INDEX>=1) THEN
               IF (PREDICTOR) TSI = T + DT
               IF (CORRECTOR) TSI = T
            ELSE
               IF (PREDICTOR) TSI = T + DT - B1%T_IGN
               IF (CORRECTOR) TSI = T      - B1%T_IGN
            ENDIF
            ! Check for cone data burning rate and compute scaled rate and time
            IF (SF%N_QDOTPP_REF > 0 .AND. N_REACTIONS>=1) THEN
               IF (PREDICTOR) THEN
                  IF (B1%T_IGN <=T) THEN
                     B1%M_DOT_G_PP_ACTUAL(N) = QDOTPP*SF%MASS_FLUX(N)
                  ELSE
                     B1%M_DOT_G_PP_ACTUAL(N) = 0._EB
                  ENDIF
               ENDIF
            ELSE
               B1%M_DOT_G_PP_ACTUAL(N) = EVALUATE_RAMP(TSI,SF%RAMP(N)%INDEX,TAU=SF%RAMP(N)%TAU)*SF%MASS_FLUX(N)
            ENDIF
            B1%M_DOT_G_PP_ADJUST(N) = B1%M_DOT_G_PP_ACTUAL(N)*SF%M_DOT_G_PP_ADJUST_FAC*B1%AREA_ADJUST
            B1%M_DOT_G_PP_ACTUAL(N) = B1%M_DOT_G_PP_ACTUAL(N)*SF%M_DOT_G_PP_ACTUAL_FAC
         ENDIF
         MFT = MFT + B1%M_DOT_G_PP_ADJUST(N)
      ENDDO SUM_MASSFLUX_LOOP

      ! Apply user-specified mass flux variation

      IF (SF%MASS_FLUX_VAR > TWO_EPSILON_EB) THEN
         ! generate pairs of standard Gaussian random variables
         CALL BOX_MULLER(RN1,RN2)
         MFT_UNIFORM = MFT
         MFT = MFT*(1._EB + RN1*SF%MASS_FLUX_VAR)
         MFT = MAX(0._EB,MIN(2._EB*MFT_UNIFORM,MFT))
         ! rescale species boundary fluxes with variation
         IF (MFT_UNIFORM>TWO_EPSILON_EB) THEN
            DO N=1,N_TRACKED_SPECIES
               B1%M_DOT_G_PP_ADJUST(N) = B1%M_DOT_G_PP_ADJUST(N) * MFT/MFT_UNIFORM
               B1%M_DOT_G_PP_ACTUAL(N) = B1%M_DOT_G_PP_ACTUAL(N) * MFT/MFT_UNIFORM
            ENDDO
         ENDIF
      ENDIF

      ! Apply water suppression coefficient (EW) at a WALL cell

      IF (PRESENT(WALL_INDEX) .OR. PRESENT(CFACE_INDEX)) THEN
         IF (B2%K_SUPPRESSION>TWO_EPSILON_EB) THEN
            B1%M_DOT_G_PP_ADJUST(1:N_TRACKED_SPECIES) = B1%M_DOT_G_PP_ADJUST(1:N_TRACKED_SPECIES)*EXP(-B2%K_SUPPRESSION)
            B1%M_DOT_G_PP_ACTUAL(1:N_TRACKED_SPECIES) = B1%M_DOT_G_PP_ACTUAL(1:N_TRACKED_SPECIES)*EXP(-B2%K_SUPPRESSION)
         ENDIF
      ENDIF

      ! Account for accumulation of ember mass (pre-transport)

      IF (CORRECTOR .AND. SF%EMBER_YIELD>0._EB .AND. (PRESENT(WALL_INDEX).OR.PRESENT(CFACE_INDEX))) THEN
         N = SF%PART_INDEX
         B2%LP_EMPUA(N) = B2%LP_EMPUA(N) + B1%M_DOT_G_PP_ACTUAL(REACTION(1)%FUEL_SMIX_INDEX)*SF%EMBER_YIELD*DT
         ! Readjust flux to conserve mass
         B1%M_DOT_G_PP_ADJUST(REACTION(1)%FUEL_SMIX_INDEX) = &
            B1%M_DOT_G_PP_ADJUST(REACTION(1)%FUEL_SMIX_INDEX)*(1._EB-SF%EMBER_YIELD)
         B1%M_DOT_G_PP_ACTUAL(REACTION(1)%FUEL_SMIX_INDEX) = &
            B1%M_DOT_G_PP_ACTUAL(REACTION(1)%FUEL_SMIX_INDEX)*(1._EB-SF%EMBER_YIELD)
      ENDIF

      ! If processing a 1-D, thermally-thick WALL cell, reduce the mass of the OBSTruction (OB%MASS) to which the WALL cell is
      ! attached. If the WALL cell is at the exterior of the current MESH, and the OBSTstruction to which it is attached
      ! lives in a neighboring MESH, store the mass to be subtracted and the index of the OBSTstruction in a 1-D array
      ! called MESHES(NM)%OMESH(NOM)%REAL_SEND_PKG8. This array will be sent to the neighboring MESH (NOM) the next time a
      ! MESH_EXCHANGE is done in main.f90.

      CONSUME_MASS: IF (PRESENT(WALL_INDEX) .AND. CORRECTOR) THEN
         OTHER_MESH_OBST_INDEX = 0
         IF (WALL_INDEX<=N_EXTERNAL_WALL_CELLS) THEN
            EWC => EXTERNAL_WALL(WALL_INDEX)
            IF (EWC%NOM>0) THEN
               IIO = EWC%IIO_MIN
               JJO = EWC%JJO_MIN
               KKO = EWC%KKO_MIN
               OTHER_MESH_OBST_INDEX = MESHES(EWC%NOM)%CELL(MESHES(EWC%NOM)%CELL_INDEX(IIO,JJO,KKO))%OBST_INDEX
            ENDIF
         ENDIF
         IF (OTHER_MESH_OBST_INDEX>0) THEN
            IF (OBST_INDEX>0) OBSTRUCTION(OBST_INDEX)%MASS = MESHES(EWC%NOM)%OBSTRUCTION(OTHER_MESH_OBST_INDEX)%MASS
            IF (MESHES(EWC%NOM)%OBSTRUCTION(OTHER_MESH_OBST_INDEX)%CONSUMABLE) THEN
               OMESH(EWC%NOM)%N_EXTERNAL_OBST = OMESH(EWC%NOM)%N_EXTERNAL_OBST + 1
               LL = 2*OMESH(EWC%NOM)%N_EXTERNAL_OBST
               OMESH(EWC%NOM)%REAL_SEND_PKG8(LL-1) = REAL(OTHER_MESH_OBST_INDEX,EB)
               OMESH(EWC%NOM)%REAL_SEND_PKG8(LL)   = &
                  (B1%M_DOT_PART_ACTUAL+SUM(B1%M_DOT_G_PP_ACTUAL(1:N_TRACKED_SPECIES)))*DT*B1%AREA
            ENDIF
         ELSE
            !$OMP CRITICAL
            IF (OBST_INDEX>0) OBSTRUCTION(OBST_INDEX)%MASS = OBSTRUCTION(OBST_INDEX)%MASS - &
               (B1%M_DOT_PART_ACTUAL+SUM(B1%M_DOT_G_PP_ACTUAL(1:N_TRACKED_SPECIES)))*DT*B1%AREA
            !$OMP END CRITICAL
         ENDIF
      ENDIF CONSUME_MASS

      ! Compute the cell face value of the species mass fraction to get the right mass flux

      IF (N_TRACKED_SPECIES==1) THEN  ! there is just the background species
         B1%RHO_F = PBAR_P(BC%KK,B1%PRESSURE_ZONE)/(RSUM0*B1%TMP_F)
         B1%ZZ_F(1) = 1._EB
         UN = MFT/B1%RHO_F
      ELSEIF (B1%LAYER_REMOVED .OR. (PRESENT(WALL_INDEX) .AND. .NOT.CELL(IC)%SOLID .AND. .NOT.CELL(IC)%EXTERIOR)) THEN
         ! this is a thin obstruction or a layer has been slated for removal
         UN = 0._EB
         B1%ZZ_F(:) = B1%ZZ_G(:)
         IF (CORRECTOR) THEN  ! calculate the mass production rate of gases in the adjacent gas cell
            RVC = RDX(BC%IIG)*RRN(BC%IIG)*RDY(BC%JJG)*RDZ(BC%KKG)
            DO NS=1,N_TRACKED_SPECIES
               IF (ABS(B1%M_DOT_G_PP_ADJUST(NS))<=TWO_EPSILON_EB) CYCLE
               MW_RATIO = SPECIES_MIXTURE(NS)%RCON/RSUM(BC%IIG,BC%JJG,BC%KKG)
               M_DOT_PPP_SINGLE = B1%M_DOT_G_PP_ADJUST(NS)*B1%AREA*RVC
               !$OMP CRITICAL
               D_SOURCE(BC%IIG,BC%JJG,BC%KKG) = D_SOURCE(BC%IIG,BC%JJG,BC%KKG) + M_DOT_PPP_SINGLE*MW_RATIO/B1%RHO_G
               M_DOT_PPP(BC%IIG,BC%JJG,BC%KKG,NS) = M_DOT_PPP(BC%IIG,BC%JJG,BC%KKG,NS) + M_DOT_PPP_SINGLE
               !$OMP END CRITICAL
            ENDDO
            ZZ_GET(1:N_TRACKED_SPECIES) = B1%ZZ_G(1:N_TRACKED_SPECIES)
            CALL GET_SPECIFIC_HEAT(ZZ_GET,CP,B1%TMP_G)
            H_G = CP*B1%TMP_G
            !$OMP CRITICAL
            D_SOURCE(BC%IIG,BC%JJG,BC%KKG) = D_SOURCE(BC%IIG,BC%JJG,BC%KKG) - B1%Q_DOT_G_PP*B1%AREA*RVC/(H_G*B1%RHO_G)
            !$OMP END CRITICAL
         ENDIF
      ELSE
         RHO_F_PREVIOUS = B1%RHO_F
         DO ITER=1,3
            UN = MFT/B1%RHO_F
            SPECIES_LOOP: DO N=1,N_TRACKED_SPECIES
               B1%RHO_D_F(N) = B1%RHO_D_F(N)*B1%RHO_F/RHO_F_PREVIOUS
               DD = 2._EB*B1%RHO_D_F(N)*B1%RDN
               DENOM = DD + UN*B1%RHO_F
               IF ( ABS(DENOM) > TWO_EPSILON_EB ) THEN
                  B1%ZZ_F(N) = ( B1%M_DOT_G_PP_ADJUST(N) + DD*B1%ZZ_G(N) ) / DENOM
               ELSE
                  B1%ZZ_F(N) = B1%ZZ_G(N)
               ENDIF
            ENDDO SPECIES_LOOP
            CALL GET_REALIZABLE_MF(B1%ZZ_F)
            ZZ_GET(1:N_TRACKED_SPECIES) = B1%ZZ_F(1:N_TRACKED_SPECIES)
            CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET,RSUM_F)
            RHO_F_PREVIOUS = B1%RHO_F
            B1%RHO_F = PBAR_P(BC%KK,B1%PRESSURE_ZONE)/(RSUM_F*B1%TMP_F)
         ENDDO
      ENDIF

      IF (PREDICTOR) B1%U_NORMAL_S = -UN
      IF (CORRECTOR) B1%U_NORMAL  = -UN

END SELECT METHOD_OF_MASS_TRANSFER

! Only set species mass fraction in the ghost cell if it is not solid

IF (PRESENT(WALL_INDEX)) THEN
   IF (WALL_INDEX<=N_EXTERNAL_WALL_CELLS .AND. .NOT.CELL(IC)%SOLID .AND. .NOT.CELL(ICG)%SOLID) &
       ZZP(BC%II,BC%JJ,BC%KK,1:N_TRACKED_SPECIES) = 2._EB*B1%ZZ_F(1:N_TRACKED_SPECIES) - B1%ZZ_G(1:N_TRACKED_SPECIES)
ENDIF

END SUBROUTINE CALCULATE_ZZ_F


!> \brief Deposit particle off-gas onto mesh
!> \param NM Mesh number
!> \param LP Pointer to Lagrangian Particle derived type variable
!> \param LPC Pointer to Lagrangian Particle Class
!> \details Deposit the particle off-gas onto the mesh

SUBROUTINE DEPOSIT_PARTICLE_MASS(NM,LP,LPC)

USE PHYSICAL_FUNCTIONS, ONLY: SURFACE_DENSITY,GET_SPECIFIC_HEAT,GET_SENSIBLE_ENTHALPY
USE OUTPUT_DATA, ONLY: M_DOT,Q_DOT
INTEGER, INTENT(IN) :: NM
REAL(EB) :: RADIUS,M_DOT_SINGLE,CP,MW_RATIO,H_G,ZZ_GET(1:N_TRACKED_SPECIES),M_GAS,LENGTH,WIDTH,H_S_B
INTEGER :: NS
TYPE(BOUNDARY_ONE_D_TYPE), POINTER :: ONE_D
TYPE(LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP
TYPE(LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC
TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC
TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1
TYPE(SURFACE_TYPE), POINTER :: SF

SF => SURFACE(LPC%SURF_INDEX)
BC => BOUNDARY_COORD(LP%BC_INDEX)
B1 => BOUNDARY_PROP1(LP%B1_INDEX)

! Get particle radius and surface area

IF (SF%PYROLYSIS_MODEL==PYROLYSIS_PREDICTED) THEN
   RADIUS = SF%INNER_RADIUS + LP%RADIUS
ELSE
   RADIUS = SF%INNER_RADIUS + SF%THICKNESS
ENDIF

IF (ABS(RADIUS)<TWO_EPSILON_EB) RETURN

IF (LPC%DRAG_LAW == SCREEN_DRAG .OR. LPC%DRAG_LAW == POROUS_DRAG) THEN
   LENGTH = LP%LENGTH
   WIDTH = LP%LENGTH
ELSE
   LENGTH = SF%LENGTH
   WIDTH = SF%WIDTH
ENDIF
SELECT CASE(SF%GEOMETRY)
   CASE(SURF_CARTESIAN)
      B1%AREA = 2._EB*LENGTH*WIDTH
   CASE(SURF_CYLINDRICAL,SURF_INNER_CYLINDRICAL)
      B1%AREA = TWOPI*RADIUS*LENGTH
   CASE(SURF_SPHERICAL)
      B1%AREA = 4._EB*PI*RADIUS**2
END SELECT

! Add evaporated particle species to gas phase and compute resulting contribution to the divergence

M_GAS = B1%RHO_G/LP%RVC
DO NS=1,N_TRACKED_SPECIES
   IF (ABS(B1%M_DOT_G_PP_ADJUST(NS))<=TWO_EPSILON_EB) CYCLE
   MW_RATIO = SPECIES_MIXTURE(NS)%RCON/RSUM(BC%IIG,BC%JJG,BC%KKG)
   M_DOT_SINGLE = LP%PWT*B1%M_DOT_G_PP_ADJUST(NS)*B1%AREA
   !$OMP CRITICAL
   D_SOURCE(BC%IIG,BC%JJG,BC%KKG) = D_SOURCE(BC%IIG,BC%JJG,BC%KKG) + M_DOT_SINGLE*(MW_RATIO/M_GAS)
   M_DOT_PPP(BC%IIG,BC%JJG,BC%KKG,NS) = M_DOT_PPP(BC%IIG,BC%JJG,BC%KKG,NS) + M_DOT_SINGLE*LP%RVC
   !$OMP END CRITICAL
   ZZ_GET = 0._EB
   ZZ_GET(NS) = 1._EB
   CALL GET_SENSIBLE_ENTHALPY(ZZ_GET,H_S_B,B1%TMP_F)
   !$OMP CRITICAL
   Q_DOT(4,NM) = Q_DOT(4,NM) + B1%M_DOT_G_PP_ADJUST(NS)*B1%AREA*H_S_B*LP%PWT    ! Q_CONV
   !$OMP END CRITICAL
ENDDO

! Calculate term in divergence that accounts for change in enthalpy between gas generated by pyrolysis and surrounding gas

ZZ_GET(1:N_TRACKED_SPECIES) = B1%ZZ_G(1:N_TRACKED_SPECIES)
CALL GET_SPECIFIC_HEAT(ZZ_GET,CP,B1%TMP_G)
H_G = CP*B1%TMP_G*M_GAS

!$OMP CRITICAL
D_SOURCE(BC%IIG,BC%JJG,BC%KKG) = D_SOURCE(BC%IIG,BC%JJG,BC%KKG) - B1%Q_DOT_G_PP*B1%AREA/H_G * LP%PWT

! Calculate contribution to divergence term due to convective heat transfer from particle

D_SOURCE(BC%IIG,BC%JJG,BC%KKG) = D_SOURCE(BC%IIG,BC%JJG,BC%KKG) - B1%Q_CON_F*B1%AREA/H_G * LP%PWT

! Add energy losses and gains to overall energy budget array

Q_DOT(8,NM) = Q_DOT(8,NM) - (B1%Q_CON_F + B1%Q_RAD_IN - B1%Q_RAD_OUT)*B1%AREA*LP%PWT      ! Q_PART
Q_DOT(3,NM) = Q_DOT(3,NM) + (B1%Q_RAD_IN-B1%Q_RAD_OUT)*B1%AREA*LP%PWT                        ! Q_RADI
!$OMP END CRITICAL

! Calculate the mass flux of fuel gas from particles

IF (CORRECTOR) THEN
   DO NS=1,N_TRACKED_SPECIES
      !$OMP CRITICAL
      M_DOT(NS,NM) = M_DOT(NS,NM) + B1%M_DOT_G_PP_ACTUAL(NS)*B1%AREA*LP%PWT
      !$OMP END CRITICAL
   ENDDO
ENDIF

! Calculate particle mass

CALC_LP_MASS: IF (SF%THERMAL_BC_INDEX==THERMALLY_THICK) THEN
   ONE_D => BOUNDARY_ONE_D(LP%OD_INDEX)
   SELECT CASE (SF%GEOMETRY)
      CASE (SURF_CARTESIAN)
         LP%MASS = B1%AREA*SF%THICKNESS*SURFACE_DENSITY(1,SF,ONE_D)
       CASE (SURF_CYLINDRICAL,SURF_INNER_CYLINDRICAL)
         LP%MASS = LENGTH*PI*(SF%INNER_RADIUS+SF%THICKNESS)**2*SURFACE_DENSITY(1,SF,ONE_D)
      CASE (SURF_SPHERICAL)
         LP%MASS = FOTHPI*(SF%INNER_RADIUS+SF%THICKNESS)**3*SURFACE_DENSITY(1,SF,ONE_D)
   END SELECT
ENDIF CALC_LP_MASS

! Update particle radius

IF (SF%PYROLYSIS_MODEL==PYROLYSIS_PREDICTED) THEN
   ONE_D => BOUNDARY_ONE_D(LP%OD_INDEX)
   LP%RADIUS = SUM(ONE_D%LAYER_THICKNESS(1:SF%N_LAYERS))
ENDIF

END SUBROUTINE DEPOSIT_PARTICLE_MASS


!> \brief Compute density, RHO_F, at non-iterpolated boundaries
!> \param BC Boundary Coordinates derived type
!> \param B1 Boundary Properties derived type
!> \param WALL_INDEX Optional WALL index
!> \param CFACE_INDEX Optional CFACE index

SUBROUTINE CALCULATE_RHO_F(BC,B1,WALL_INDEX,CFACE_INDEX)

USE PHYSICAL_FUNCTIONS, ONLY : GET_SPECIFIC_GAS_CONSTANT
INTEGER, INTENT(IN), OPTIONAL :: WALL_INDEX,CFACE_INDEX
TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1
TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC
REAL(EB) :: ZZ_GET(1:N_TRACKED_SPECIES),RSUM_F,UN_P
INTEGER :: BOUNDARY_TYPE

IF (PRESENT(WALL_INDEX)) THEN
   BOUNDARY_TYPE = WALL(WALL_INDEX)%BOUNDARY_TYPE
ELSEIF (PRESENT(CFACE_INDEX)) THEN
   BOUNDARY_TYPE = CFACE(CFACE_INDEX)%BOUNDARY_TYPE
ENDIF

ZZ_GET(1:N_TRACKED_SPECIES) = MAX(0._EB,B1%ZZ_F(1:N_TRACKED_SPECIES))
CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET,RSUM_F)
B1%RHO_F = PBAR_P(BC%KK,B1%PRESSURE_ZONE)/(RSUM_F*B1%TMP_F)

! If the boundary is solid and gas is being drawn in, set surface variables to equal the adjacent gas phase variables

IF (BOUNDARY_TYPE==SOLID_BOUNDARY) THEN
   IF (PREDICTOR) THEN
      UN_P = B1%U_NORMAL_S
   ELSE
      UN_P = B1%U_NORMAL
   ENDIF
   IF (UN_P>0._EB) THEN
      B1%ZZ_F(1:N_TRACKED_SPECIES) = B1%ZZ_G(1:N_TRACKED_SPECIES)
      B1%RHO_F = B1%RHO_G
   ENDIF
ENDIF

IF (PRESENT(WALL_INDEX)) THEN
   IF (WALL_INDEX<=N_EXTERNAL_WALL_CELLS .AND. BOUNDARY_TYPE/=OPEN_BOUNDARY) RHOP(BC%II,BC%JJ,BC%KK) = 2._EB*B1%RHO_F - B1%RHO_G
ENDIF

END SUBROUTINE CALCULATE_RHO_F


!> \brief Calculate aerosol deposition onto a solid surface
!> \param DT Current time step
!> \param BC Pointer to Boundary Coordinate derived type variable
!> \param B1 Pointer to Boundary Property derived type variable
!> \param B2 Pointer to Boundary Property derived type variable
!> \param WALL_INDEX Optional WALL index
!> \param CFACE_INDEX Optional CFACE index

SUBROUTINE CALC_DEPOSITION(DT,BC,B1,B2,WALL_INDEX,CFACE_INDEX)

USE PHYSICAL_FUNCTIONS, ONLY: GET_VISCOSITY,GET_CONDUCTIVITY,CUNNINGHAM
USE GLOBAL_CONSTANTS, ONLY: K_BOLTZMANN,GRAVITATIONAL_DEPOSITION,TURBULENT_DEPOSITION,THERMOPHORETIC_DEPOSITION,GVEC
REAL(EB), INTENT(IN) :: DT
INTEGER, INTENT(IN), OPTIONAL :: WALL_INDEX,CFACE_INDEX
INTEGER:: N
TYPE(SPECIES_MIXTURE_TYPE), POINTER :: SM
TYPE(SPECIES_TYPE), POINTER :: SS
REAL(EB), PARAMETER :: CS=1.17_EB,CT=2.2_EB,CM=1.146_EB
REAL(EB), PARAMETER :: CM3=3._EB*CM,CS2=CS*2._EB,CT2=2._EB*CT
REAL(EB), PARAMETER :: ZZ_MIN_DEP=1.E-14_EB
REAL(EB) :: U_THERM,U_TURB,MU_G,Y_AEROSOL,ZZ_GET(1:N_TRACKED_SPECIES),YDEP,K_G,TMP_FILM,ALPHA,DTMPDX,&
            TAU_PLUS,U_GRAV,D_SOLID,MW_RATIO,KN,KN_FAC,RSUM_G,ZZ_G(1:N_TRACKED_SPECIES),TAU_PLUS_C
TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC
TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1
TYPE(BOUNDARY_PROP2_TYPE), POINTER :: B2

IF (PRESENT(WALL_INDEX)) THEN
   RSUM_G = RSUM(BC%IIG,BC%JJG,BC%KKG)
   ZZ_G(1:N_TRACKED_SPECIES) = ZZ(BC%IIG,BC%JJG,BC%KKG,1:N_TRACKED_SPECIES)
ELSEIF (PRESENT(CFACE_INDEX)) THEN
   RSUM_G= CFACE(CFACE_INDEX)%RSUM_G
   ZZ_G(1:N_TRACKED_SPECIES) = B1%ZZ_G(1:N_TRACKED_SPECIES)
ENDIF

SMIX_LOOP: DO N=1,N_TRACKED_SPECIES

   ZZ_GET(1:N_TRACKED_SPECIES) = MAX(0._EB,ZZ_G(1:N_TRACKED_SPECIES))
   IF (ZZ_GET(N) < ZZ_MIN_DEP) CYCLE SMIX_LOOP
   SM => SPECIES_MIXTURE(N)
   IF (.NOT.SM%DEPOSITING) CYCLE SMIX_LOOP
   SS => SPECIES(SM%SINGLE_SPEC_INDEX)
   MW_RATIO = SPECIES_MIXTURE(N)%RCON/RSUM_G
   TMP_FILM = 0.5_EB*(B1%TMP_G+B1%TMP_F)
   CALL GET_VISCOSITY(ZZ_GET,MU_G,TMP_FILM)
   CALL GET_CONDUCTIVITY(ZZ_GET,K_G,TMP_FILM)
   KN_FAC = MU_G*SQRT(0.5_EB*PI/(PBAR(BC%KKG,PRESSURE_ZONE(BC%IIG,BC%JJG,BC%KKG))*B1%RHO_G))
   ALPHA = K_G/SM%CONDUCTIVITY_SOLID
   DTMPDX = B1%HEAT_TRANS_COEF*(B1%TMP_G-B1%TMP_F)/K_G
   U_THERM = 0._EB
   U_TURB = 0._EB
   U_GRAV = 0._EB

   IF (THERMOPHORETIC_DEPOSITION) THEN
      KN = KN_FAC/SM%THERMOPHORETIC_DIAMETER
      U_THERM = CS2*(ALPHA+CT*KN)*CUNNINGHAM(KN)/((1._EB+CM3*KN)*(1+2*ALPHA+CT2*KN)) * MU_G/(B1%TMP_G*B1%RHO_G)*DTMPDX
   ENDIF

   IF (GRAVITATIONAL_DEPOSITION) THEN
      KN = KN_FAC/SM%MEAN_DIAMETER
      U_GRAV = - DOT_PRODUCT(GVEC,BC%NVEC)*CUNNINGHAM(KN)*SM%MEAN_DIAMETER**2*SM%DENSITY_SOLID/(18._EB*MU_G)
      U_GRAV = MAX(0._EB,U_GRAV)  ! Prevent negative settling velocity at downward facing surfaces
   ENDIF

   IF (TURBULENT_DEPOSITION) THEN
      KN = KN_FAC/SM%MEAN_DIAMETER
      TAU_PLUS_C = SM%DENSITY_SOLID*SM%MEAN_DIAMETER**2/18._EB
      TAU_PLUS = TAU_PLUS_C/MU_G**2*B2%U_TAU**2*B1%RHO_G
      IF (TAU_PLUS < 0.2_EB) THEN ! Diffusion regime
         D_SOLID = K_BOLTZMANN*B1%TMP_G*CUNNINGHAM(KN)/(3._EB*PI*MU_G*SM%MEAN_DIAMETER)
         U_TURB = B2%U_TAU * 0.086_EB*(MU_G/B1%RHO_G/D_SOLID)**(-0.7_EB)
      ELSEIF (TAU_PLUS >= 0.2_EB .AND. TAU_PLUS < 22.9_EB) THEN ! Diffusion-impaction regime
         U_TURB = B2%U_TAU * 3.5E-4_EB * TAU_PLUS**2
      ELSE ! Inertia regime
         U_TURB = B2%U_TAU * 0.17_EB
      ENDIF
   ENDIF

   B2%V_DEP = MAX(0._EB,U_THERM+U_TURB+U_GRAV+B1%U_NORMAL)
   IF (B2%V_DEP <= TWO_EPSILON_EB) CYCLE SMIX_LOOP
   ZZ_GET = ZZ_GET * B1%RHO_G
   Y_AEROSOL = ZZ_GET(N)
   YDEP = Y_AEROSOL*MIN(1._EB,(B2%V_DEP)*DT*B1%RDN)
   ZZ_GET(N) = Y_AEROSOL - YDEP
   IF (SM%AWM_INDEX > 0) B1%AWM_AEROSOL(SM%AWM_INDEX)= B1%AWM_AEROSOL(SM%AWM_INDEX)+YDEP/B1%RDN
   IF (SS%AWM_INDEX > 0) B1%AWM_AEROSOL(SS%AWM_INDEX)= B1%AWM_AEROSOL(SS%AWM_INDEX)+YDEP/B1%RDN

   !$OMP CRITICAL
   D_SOURCE(BC%IIG,BC%JJG,BC%KKG) = D_SOURCE(BC%IIG,BC%JJG,BC%KKG) - MW_RATIO*YDEP / B1%RHO_G / DT
   M_DOT_PPP(BC%IIG,BC%JJG,BC%KKG,N) = M_DOT_PPP(BC%IIG,BC%JJG,BC%KKG,N) - YDEP / DT
   !$OMP END CRITICAL

ENDDO SMIX_LOOP

END SUBROUTINE CALC_DEPOSITION


!> \brief Calculate HVAC boundary conditions
!> \param BC Pointer to Boundary Coordinate derived type variable
!> \param B1 Pointer to Boundary Property derived type variable
!> \param SF Pointer to Surface Propertes derived type variable

SUBROUTINE CALC_HVAC_BC(BC,B1,SF)

USE HVAC_ROUTINES, ONLY : NODE_AREA_EX,NODE_TMP_EX,DUCT_MF,NODE_ZZ_EX
USE PHYSICAL_FUNCTIONS, ONLY : GET_SPECIFIC_GAS_CONSTANT,GET_ENTHALPY
REAL(EB) :: ZZ_GET(1:N_TRACKED_SPECIES),UN,MFT,RSUM_F,H_D,H_G
INTEGER  :: COUNTER,DU
TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1
TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC
TYPE(SURFACE_TYPE), POINTER :: SF

COUNTER = 0

! Compute R*Sum(Y_i/W_i) at the wall

DU=DUCTNODE(B1%NODE_INDEX)%DUCT_INDEX(1)
MFT = -DUCTNODE(B1%NODE_INDEX)%DIR(1)*DUCT_MF(DU)/NODE_AREA_EX(B1%NODE_INDEX)
IF (.NOT. ANY(SF%LEAK_PATH>0)) THEN
   IF (DUCTNODE(B1%NODE_INDEX)%DIR(1)*DUCT_MF(DU) > 0._EB) THEN
      IF (SF%THERMAL_BC_INDEX==HVAC_BOUNDARY) THEN
         B1%TMP_F = NODE_TMP_EX(B1%NODE_INDEX)
         B1%HEAT_TRANS_COEF = 0._EB
         B1%Q_CON_F = 0._EB
      ELSE
         IF (DUCT(DU)%LEAK_ENTHALPY) THEN
            ZZ_GET(1:N_TRACKED_SPECIES) = NODE_ZZ_EX(B1%NODE_INDEX,1:N_TRACKED_SPECIES)
            CALL GET_ENTHALPY(ZZ_GET,H_G,B1%TMP_F)
            CALL GET_ENTHALPY(ZZ_GET,H_D,NODE_TMP_EX(B1%NODE_INDEX))
            B1%Q_LEAK = -MFT*(H_D-H_G)*B1%RDN
         ENDIF
      ENDIF
   ELSE
      IF (SF%THERMAL_BC_INDEX==HVAC_BOUNDARY) THEN
         B1%TMP_F = B1%TMP_G
         B1%HEAT_TRANS_COEF = 0._EB
         B1%Q_CON_F = 0._EB
      ENDIF
   ENDIF
ENDIF

IF (MFT >= 0._EB) THEN
   ZZ_GET(1:N_TRACKED_SPECIES) = B1%ZZ_G(1:N_TRACKED_SPECIES)
   CALL GET_SPECIFIC_GAS_CONSTANT(ZZ_GET,RSUM_F)
   B1%RHO_F = PBAR_P(BC%KK,B1%PRESSURE_ZONE)/(RSUM_F*B1%TMP_G)
   UN = MFT/B1%RHO_F
   IF (PREDICTOR) B1%U_NORMAL_S = UN
   IF (CORRECTOR) B1%U_NORMAL  = UN
ELSE
   B1%M_DOT_G_PP_ADJUST(1:N_TRACKED_SPECIES) = -NODE_ZZ_EX(B1%NODE_INDEX,1:N_TRACKED_SPECIES)*MFT
ENDIF

END SUBROUTINE CALC_HVAC_BC


!> \brief Update temperature and material components for a single boundary cell
!>
!> \param NM Mesh number
!> \param T Current time (s)
!> \param DT_BC Time step (s) used for solid phase updates
!> \param PARTICLE_INDEX Index of a Lagrangian particle
!> \param WALL_INDEX Index of a Cartesian WALL cell
!> \param CFACE_INDEX Index of an immersed boundary CFACE

SUBROUTINE SOLID_HEAT_TRANSFER(NM,T,DT_BC,PARTICLE_INDEX,WALL_INDEX,CFACE_INDEX,THIN_WALL_INDEX)

USE GEOMETRY_FUNCTIONS, ONLY: GET_WALL_NODE_WEIGHTS,GET_N_LAYER_CELLS,GET_WALL_NODE_COORDINATES,INTERPOLATE_WALL_ARRAY,&
                              GET_INTERPOLATION_WEIGHTS
USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP,INTERPOLATE1D_UNIFORM
USE COMP_FUNCTIONS, ONLY: SHUTDOWN
USE PHYSICAL_FUNCTIONS, ONLY: GET_SPECIFIC_GAS_CONSTANT,GET_SENSIBLE_ENTHALPY
REAL(EB), INTENT(IN) :: DT_BC,T
INTEGER, INTENT(IN) :: NM
INTEGER, INTENT(IN), OPTIONAL:: WALL_INDEX,PARTICLE_INDEX,CFACE_INDEX,THIN_WALL_INDEX
REAL(EB) :: RDT,DTMP,QDXKF,QDXKB,RR,RFACF,RFACB,RFACF2,RFACB2,Q_RAD_IN_B,RFLUX_UP,RFLUX_DOWN,E_WALLB, &
            VOLSUM,KAPSUM,DXF,DXB,HTCF,HTCB,Q_RAD_OUT,Q_RAD_OUT_OLD,Q_CON_F,Q_CON_B,&
            Q_LIQUID_F,Q_LIQUID_B,LAYER_DIVIDE,TMP_GAS_BACK,GEOM_FACTOR,DT_BC_SUB_OLD,&
            DEL_DOT_Q_SC,Q_DOT_G_PP,Q_DOT_G_PP_NET,Q_DOT_O2_PP,Q_DOT_O2_PP_NET,R_SURF,U_SURF,V_SURF,W_SURF,T_BC_SUB,DT_BC_SUB,&
            Q_NET_F,Q_NET_B,TMP_RATIO,KODXF,KODXB,H_S,T_NODE,C_S,H_NODE,VOL,&
            RADIUS,HTC_LIMIT,CP1,CP2,DENOM,SF_HTC_F,SF_HTC_B,THICKNESS,DT_FO,DDSUM,NODE_RDT(NWP_MAX)
REAL(EB), DIMENSION(N_TRACKED_SPECIES) :: M_DOT_G_PP_ADJUST,M_DOT_G_PP_ADJUST_NET,M_DOT_G_PP_ACTUAL,M_DOT_G_PP_ACTUAL_NET
REAL(EB), DIMENSION(MAX_MATERIALS) :: M_DOT_S_PP,M_DOT_S_PP_NET,T_BOIL_EFF
REAL(EB), DIMENSION(MAX_LPC) :: Q_DOT_PART_S,M_DOT_PART_S
REAL(EB), DIMENSION(NWP_MAX) :: TMP_S,RHO_H_S
REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: RHO_DOT,INT_WGT
REAL(EB), DIMENSION(MAX_LAYERS) :: DX_MIN
REAL(EB), DIMENSION(MAX_LAYERS,MAX_MATERIALS) :: RHO_ADJUSTED
REAL(EB), DIMENSION(NWP_MAX) :: AAS,BBS,CCS,DDS,DDT,Q_S,Q_IR,Q_ADD,TWO_DX_KAPPA_S,DX_S,MF_FRAC,REGRID_FACTOR
REAL(EB), DIMENSION(0:NWP_MAX+1) :: RHO_S,DELTA_TMP,RDX_S
REAL(EB), DIMENSION(0:NWP_MAX) :: X_S_NEW,RDXN_S,R_S,R_S_NEW,DX_WGT_S
INTEGER, DIMENSION(0:NWP_MAX+1) :: LAYER_INDEX
INTEGER, DIMENSION(MAX_LAYERS) :: N_LAYER_CELLS_NEW
INTEGER  :: NWP_NEW,I_GRAD,IZERO,SURF_INDEX,BACKING,NWP,I,NL,N,OBST_INDEX,&
            N_CELLS,ITMP,ITER,BACK_MESH,BACK_INDEX,BACK_WALL_INDEX
CHARACTER(MESSAGE_LENGTH) :: MESSAGE
LOGICAL :: ISOLATED_THIN_WALL,ISOLATED_THIN_WALL_BACK,E_FOUND,REMESH_LAYER(MAX_LAYERS),REMESH_CHECK,&
           CELL_ZERO(MAX_LAYERS),TMP_CHECK(MAX_LAYERS)
TYPE(WALL_TYPE), POINTER :: WC,WC_BACK
TYPE(THIN_WALL_TYPE), POINTER :: TW,TW_BACK
TYPE(CFACE_TYPE), POINTER :: CFA,CFA_BACK
TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1,B1M,B1P,B1_BACK
TYPE(BOUNDARY_PROP2_TYPE), POINTER :: B2,B2_BACK
TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC,BC_BACK
TYPE(BOUNDARY_ONE_D_TYPE), POINTER :: ONE_D
TYPE(LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP
TYPE(SURFACE_TYPE), POINTER :: SF,SF_BACK
TYPE(MATERIAL_TYPE), POINTER :: ML

! Copy commonly used derived type variables into local variables.

R_SURF=0._EB
U_SURF=0._EB
V_SURF=0._EB
W_SURF=0._EB

ISOLATED_THIN_WALL = .FALSE.
ISOLATED_THIN_WALL_BACK = .FALSE.

UNPACK_WALL_PARTICLE: IF (PRESENT(WALL_INDEX)) THEN

   WC => WALL(WALL_INDEX)
   SURF_INDEX = WC%SURF_INDEX
   SF => SURFACE(SURF_INDEX)
   ONE_D => BOUNDARY_ONE_D(WC%OD_INDEX)
   B1 => BOUNDARY_PROP1(WC%B1_INDEX)
   B2 => BOUNDARY_PROP2(WC%B2_INDEX)
   BC => BOUNDARY_COORD(WC%BC_INDEX)
   IF (SF%HT_DIM>1 .AND. ABS(BC%IOR)/=HT_3D_SWEEP_DIRECTION) RETURN
   OBST_INDEX      = WC%OBST_INDEX
   BACK_INDEX      = ONE_D%BACK_INDEX
   BACK_WALL_INDEX = BACK_INDEX
   BACK_MESH       = ONE_D%BACK_MESH
   IF (BACK_INDEX>0) THEN
      WC_BACK => MESHES(BACK_MESH)%WALL(BACK_INDEX)
      B1_BACK => MESHES(BACK_MESH)%BOUNDARY_PROP1(WC_BACK%B1_INDEX)
      B2_BACK => MESHES(BACK_MESH)%BOUNDARY_PROP2(WC_BACK%B2_INDEX)
      BC_BACK => MESHES(BACK_MESH)%BOUNDARY_COORD(WC_BACK%BC_INDEX)
      SF_BACK => SURFACE(WC_BACK%SURF_INDEX)
   ENDIF
   BACKING = SF%BACKING
   IF (BACK_INDEX==0 .AND. SF%BACKING==EXPOSED) BACKING = VOID

   ! Energy flux due to evaporation or condensation

   IF (LIQUID_DROPLETS) THEN
      Q_LIQUID_F = -SUM(B2%LP_CPUA(:)) + B1%Q_CONDENSE
   ELSE
      Q_LIQUID_F = B1%Q_CONDENSE
   ENDIF

ELSEIF (PRESENT(THIN_WALL_INDEX)) THEN

   TW => THIN_WALL(THIN_WALL_INDEX)
   SURF_INDEX = TW%SURF_INDEX
   SF => SURFACE(SURF_INDEX)
   ONE_D => BOUNDARY_ONE_D(TW%OD_INDEX)
   B1 => BOUNDARY_PROP1(TW%B1_INDEX)
   BC => BOUNDARY_COORD(TW%BC_INDEX)
   IF (SF%HT_DIM>1 .AND. ABS(BC%IOR)/=HT_3D_SWEEP_DIRECTION) RETURN
   OBST_INDEX = TW%OBST_INDEX
   BACK_INDEX = ONE_D%BACK_INDEX
   BACK_MESH  = ONE_D%BACK_MESH
   IF (MESHES(BACK_MESH)%THIN_WALL(BACK_INDEX)%WALL_INDEX_M>0) THEN
      BACK_WALL_INDEX = MESHES(BACK_MESH)%THIN_WALL(BACK_INDEX)%WALL_INDEX_M
   ELSEIF (MESHES(BACK_MESH)%THIN_WALL(BACK_INDEX)%WALL_INDEX_P>0) THEN
      BACK_WALL_INDEX = MESHES(BACK_MESH)%THIN_WALL(BACK_INDEX)%WALL_INDEX_P
   ELSE
      BACK_WALL_INDEX = 0
      ISOLATED_THIN_WALL_BACK = .TRUE.
   ENDIF
   IF (BACK_INDEX>0) THEN
      TW_BACK => MESHES(BACK_MESH)%THIN_WALL(BACK_INDEX)
      B1_BACK => MESHES(BACK_MESH)%BOUNDARY_PROP1(TW_BACK%B1_INDEX)
      BC_BACK => MESHES(BACK_MESH)%BOUNDARY_COORD(TW_BACK%BC_INDEX)
      SF_BACK => SURFACE(TW_BACK%SURF_INDEX)
   ENDIF
   IF (BACK_WALL_INDEX>0) B1_BACK => MESHES(BACK_MESH)%BOUNDARY_PROP1(MESHES(BACK_MESH)%WALL(BACK_WALL_INDEX)%B1_INDEX)
   BACKING = SF%BACKING
   Q_LIQUID_F = 0._EB
   Q_LIQUID_B = 0._EB
   IF (TW%WALL_INDEX_M>0 .OR. TW%WALL_INDEX_P>0) THEN
      IF (TW%WALL_INDEX_M>0) THEN
         B1M => BOUNDARY_PROP1(WALL(TW%WALL_INDEX_M)%B1_INDEX)
      ELSE
         B1M => BOUNDARY_PROP1(WALL(TW%WALL_INDEX_P)%B1_INDEX)
      ENDIF
      IF (TW%WALL_INDEX_P>0) THEN
         B1P => BOUNDARY_PROP1(WALL(TW%WALL_INDEX_P)%B1_INDEX)
      ELSE
         B1P => BOUNDARY_PROP1(WALL(TW%WALL_INDEX_M)%B1_INDEX)
      ENDIF
      B1%TMP_G = 0.5_EB*(B1M%TMP_G+B1P%TMP_G)
      B1%Q_RAD_IN = 0.5_EB*(B1M%Q_RAD_IN+B1P%Q_RAD_IN)
   ELSE
      B1%TMP_G = TMP(BC%IIG,BC%JJG,BC%KKG)
      ISOLATED_THIN_WALL = .TRUE.
   ENDIF
   ! Special case where the gas temperature is fixed by the user

   IF (SF%TMP_GAS_FRONT > 0._EB) B1%TMP_G = TMPA + EVALUATE_RAMP(T-T_BEGIN,SF%RAMP(TIME_TGF)%INDEX)*(SF%TMP_GAS_FRONT-TMPA)

ELSEIF (PRESENT(CFACE_INDEX)) THEN UNPACK_WALL_PARTICLE

   CFA => CFACE(CFACE_INDEX)
   SURF_INDEX = CFA%SURF_INDEX
   SF => SURFACE(SURF_INDEX)
   ONE_D => BOUNDARY_ONE_D(CFA%OD_INDEX)
   B1 => BOUNDARY_PROP1(CFA%B1_INDEX)
   B2 => BOUNDARY_PROP2(CFA%B2_INDEX)
   BC => BOUNDARY_COORD(CFA%BC_INDEX)
   OBST_INDEX = 0
   BACK_INDEX      = ONE_D%BACK_INDEX
   BACK_WALL_INDEX = BACK_INDEX
   BACK_MESH       = ONE_D%BACK_MESH
   IF (LIQUID_DROPLETS) THEN
      Q_LIQUID_F = -SUM(B2%LP_CPUA(:)) + B1%Q_CONDENSE
   ELSE
      Q_LIQUID_F = B1%Q_CONDENSE
   ENDIF
   IF (BACK_INDEX>0 .AND. BACK_MESH==NM) THEN
      CFA_BACK => CFACE(BACK_INDEX)
      B1_BACK => BOUNDARY_PROP1(CFA_BACK%B1_INDEX)
      B2_BACK => BOUNDARY_PROP2(CFA_BACK%B2_INDEX)
      BC_BACK => BOUNDARY_COORD(CFA_BACK%BC_INDEX)
      SF_BACK => SURFACE(CFA_BACK%SURF_INDEX)
   ENDIF
   BACKING = SF%BACKING
   IF (ONE_D%BACK_INDEX==0 .AND. SF%BACKING==EXPOSED) BACKING = VOID

ELSEIF (PRESENT(PARTICLE_INDEX)) THEN UNPACK_WALL_PARTICLE

   LP => LAGRANGIAN_PARTICLE(PARTICLE_INDEX)
   ONE_D => BOUNDARY_ONE_D(LP%OD_INDEX)
   B1 => BOUNDARY_PROP1(LP%B1_INDEX)
   IF (OXPYRO_MODEL) B2 => BOUNDARY_PROP2(LP%B2_INDEX)
   BC => BOUNDARY_COORD(LP%BC_INDEX)
   SURF_INDEX = LAGRANGIAN_PARTICLE_CLASS(LP%CLASS_INDEX)%SURF_INDEX
   SF => SURFACE(SURF_INDEX)
   OBST_INDEX = 0
   BACK_INDEX      = 0
   BACK_WALL_INDEX = BACK_INDEX
   BACK_MESH       = NM
   IF (LAGRANGIAN_PARTICLE_CLASS(LP%CLASS_INDEX)%MASSLESS_TARGET .AND. SF%GEOMETRY==SURF_CARTESIAN) THEN
      BACKING = VOID
   ELSE
      BACKING = INSULATED
   ENDIF
   Q_LIQUID_F = B1%Q_CONDENSE
   R_SURF=SUM(ONE_D%LAYER_THICKNESS(1:ONE_D%N_LAYERS))
   U_SURF=LP%U
   V_SURF=LP%V
   W_SURF=LP%W

ENDIF UNPACK_WALL_PARTICLE

! Assume that no layers have been slated for removal

B1%LAYER_REMOVED = .FALSE.

! Evaluate any heat transfer coefficient ramp

IF (SF%H_FIXED >= 0._EB) THEN
   SF_HTC_F = SF%H_FIXED
   IF (SF%RAMP_H_FIXED_INDEX > 0) SF_HTC_F = SF_HTC_F * EVALUATE_RAMP(T-T_BEGIN,SF%RAMP_H_FIXED_INDEX)
ELSE
   SF_HTC_F = -1._EB
ENDIF

IF (SF%H_FIXED_B >= 0._EB) THEN
   SF_HTC_B = SF%H_FIXED_B
   IF (SF%RAMP_H_FIXED_B_INDEX > 0) SF_HTC_B = SF_HTC_B * EVALUATE_RAMP(T-T_BEGIN,SF%RAMP_H_FIXED_B_INDEX)
ELSE
   SF_HTC_B = -1._EB
ENDIF

! If the fuel has burned away, return

IF (B1%BURNAWAY) THEN
   B1%M_DOT_G_PP_ADJUST(1:N_TRACKED_SPECIES) = 0._EB
   B1%M_DOT_G_PP_ACTUAL(1:N_TRACKED_SPECIES) = 0._EB
   ONE_D%M_DOT_S_PP(1:ONE_D%N_MATL) = 0._EB
   RETURN
ENDIF

! Exponents for cylindrical or spherical coordinates

SELECT CASE(SF%GEOMETRY)
   CASE(SURF_CARTESIAN)                          ; I_GRAD = 1
   CASE(SURF_CYLINDRICAL,SURF_INNER_CYLINDRICAL) ; I_GRAD = 2
   CASE(SURF_SPHERICAL)                          ; I_GRAD = 3
END SELECT

! Create array of material densities that are adjusted in special cases

DO NL=1,ONE_D%N_LAYERS
   DO N=1,ONE_D%N_MATL
      RHO_ADJUSTED(NL,N) = MATERIAL(ONE_D%MATL_INDEX(N))%RHO_S*SF%DENSITY_ADJUST_FACTOR(NL,N)
   ENDDO
ENDDO

! Set mass and energy fluxes to zero prior to time sub-iteration

B1%Q_CON_F = 0._EB
IF (ONE_D%INTERNAL_RADIATION) Q_RAD_OUT_OLD = B1%Q_RAD_OUT
B1%Q_RAD_OUT = 0._EB

IF (ONE_D%PYROLYSIS_MODEL==PYROLYSIS_PREDICTED) THEN

   IF (ALLOCATED(RHO_DOT)) DEALLOCATE(RHO_DOT)
   ALLOCATE(RHO_DOT(ONE_D%N_MATL,NWP_MAX))

   M_DOT_G_PP_ADJUST_NET          = 0._EB
   M_DOT_G_PP_ACTUAL_NET          = 0._EB
   M_DOT_S_PP_NET(1:ONE_D%N_MATL) = 0._EB
   Q_DOT_G_PP_NET                 = 0._EB
   Q_DOT_O2_PP_NET                = 0._EB
   B1%M_DOT_PART_ACTUAL           = 0._EB
ENDIF

! Start time iterations here

T_BC_SUB  = 0._EB
DT_BC_SUB = DT_BC
B1%N_SUBSTEPS = 1
ONE_D%DELTA_TMP = 0._EB

! Compute initial thermal properties for Fo number calculation
! CHECK_FO=F by default
! CHECK_FO=T enforces near explicit time step accuracy (Forward-Euler would require 1/2 DT_FO)

CHECK_FO_IF: IF (CHECK_FO) THEN
   ONE_D%K_S = 0._EB
   RHO_S   = 0._EB
   ONE_D%RHO_C_S = 0._EB
   POINT_LOOP0: DO I=1,NWP
      VOLSUM = 0._EB
      ITMP = MIN(I_MAX_TEMP-1,INT(ONE_D%TMP(I)))
      MATERIAL_LOOP0: DO N=1,ONE_D%N_MATL
         IF (ONE_D%MATL_COMP(N)%RHO(I)<=TWO_EPSILON_EB) CYCLE MATERIAL_LOOP0
         ML  => MATERIAL(ONE_D%MATL_INDEX(N))
         VOLSUM = VOLSUM + ONE_D%MATL_COMP(N)%RHO(I)/RHO_ADJUSTED(LAYER_INDEX(I),N)
         ONE_D%K_S(I) = ONE_D%K_S(I) + ONE_D%MATL_COMP(N)%RHO(I)*ML%K_S(ITMP)/RHO_ADJUSTED(LAYER_INDEX(I),N)
         ONE_D%RHO_C_S(I) = ONE_D%RHO_C_S(I) + ONE_D%MATL_COMP(N)%RHO(I)*ML%C_S(ITMP)
         RHO_S(I) = RHO_S(I) + ONE_D%MATL_COMP(N)%RHO(I)
      ENDDO MATERIAL_LOOP0
      IF (SF%PACKING_RATIO(LAYER_INDEX(I))>0._EB) ONE_D%K_S(I) = ONE_D%K_S(I)*SF%PACKING_RATIO(LAYER_INDEX(I))

      IF (VOLSUM > 0._EB) THEN
         ONE_D%K_S(I) = ONE_D%K_S(I)/VOLSUM
      ENDIF
      IF (ONE_D%K_S(I)<=TWO_EPSILON_EB)      ONE_D%K_S(I)      = 10000._EB
      IF (ONE_D%RHO_C_S(I)<=TWO_EPSILON_EB)  ONE_D%RHO_C_S(I)  = 0.001_EB
   ENDDO POINT_LOOP0
ENDIF CHECK_FO_IF

SUB_TIMESTEP_LOOP: DO

   ! Compute grid for reacting nodes

   LAYER_DIVIDE = SF%LAYER_DIVIDE

   COMPUTE_GRID: IF (ONE_D%PYROLYSIS_MODEL==PYROLYSIS_PREDICTED .OR. SF%HT_DIM>1 .OR. SF%VARIABLE_THICKNESS) THEN
      NWP = SUM(ONE_D%N_LAYER_CELLS(1:ONE_D%N_LAYERS))
      CALL GET_WALL_NODE_WEIGHTS(NWP,ONE_D%N_LAYERS,ONE_D%N_LAYER_CELLS(1:ONE_D%N_LAYERS),ONE_D%LAYER_THICKNESS,SF%GEOMETRY, &
         ONE_D%X(0:NWP),LAYER_DIVIDE,DX_S(1:NWP),RDX_S(0:NWP+1),RDXN_S(0:NWP),DX_WGT_S(0:NWP),DXF,DXB,&
         LAYER_INDEX(0:NWP+1),MF_FRAC(1:NWP),SF%INNER_RADIUS)
   ELSE COMPUTE_GRID
      NWP                  = SF%N_CELLS_INI
      DXF                  = SF%DXF
      DXB                  = SF%DXB
      DX_S(1:NWP)          = SF%DX(1:NWP)
      RDX_S(0:NWP+1)       = SF%RDX(0:NWP+1)
      RDXN_S(0:NWP)        = SF%RDXN(0:NWP)
      DX_WGT_S(0:NWP)      = SF%DX_WGT(0:NWP)
      LAYER_INDEX(0:NWP+1) = SF%LAYER_INDEX(0:NWP+1)
      MF_FRAC(1:NWP)       = SF%MF_FRAC(1:NWP)
   ENDIF COMPUTE_GRID

   ! Calculate minimum DT based on Fourier number, Fo

   IF (CHECK_FO) THEN
      DT_FO = MINVAL( DX_S(1:NWP)**2 * ONE_D%RHO_C_S(1:NWP) / ONE_D%K_S(1:NWP) )
   ELSE
      DT_FO = HUGE(1._EB)
   ENDIF

   ! Compute convective heat flux at the surface

   DTMP = B1%TMP_G - B1%TMP_F
   IF (PRESENT(WALL_INDEX)) THEN
      HTCF = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC_F,SF,WALL_INDEX_IN=WALL_INDEX)
   ELSEIF (PRESENT(THIN_WALL_INDEX)) THEN
      IF (TW%WALL_INDEX_M>0) THEN
         HTCF = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC_F,SF,WALL_INDEX_IN=TW%WALL_INDEX_M)
      ELSEIF (TW%WALL_INDEX_P>0) THEN
         HTCF = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC_F,SF,WALL_INDEX_IN=TW%WALL_INDEX_P)
      ELSE
         HTCF = 0._EB
      ENDIF
   ELSEIF (PRESENT(CFACE_INDEX)) THEN
      HTCF = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC_F,SF,CFACE_INDEX_IN=CFACE_INDEX)
   ELSEIF (PRESENT(PARTICLE_INDEX)) THEN
      RADIUS = SF%INNER_RADIUS + SUM(ONE_D%LAYER_THICKNESS(1:SF%N_LAYERS))
      SELECT CASE(SF%GEOMETRY)
         CASE (SURF_CARTESIAN)                          ; HTC_LIMIT = 0.5_EB*RADIUS*ONE_D%RHO_C_S(1)/(      DT_BC_SUB)
         CASE (SURF_CYLINDRICAL,SURF_INNER_CYLINDRICAL) ; HTC_LIMIT = 0.5_EB*RADIUS*ONE_D%RHO_C_S(1)/(2._EB*DT_BC_SUB)
         CASE (SURF_SPHERICAL)                          ; HTC_LIMIT = 0.5_EB*RADIUS*ONE_D%RHO_C_S(1)/(3._EB*DT_BC_SUB)
      END SELECT
      HTCF = MIN(HTC_LIMIT , HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC_F,SF,PARTICLE_INDEX_IN=PARTICLE_INDEX))
   ENDIF
   Q_CON_F = HTCF*DTMP

   ! Compute back side emissivity

   IF (SF%EMISSIVITY_BACK_SPECIFIED) THEN
      E_WALLB = SF%EMISSIVITY_BACK
   ELSEIF (BACKING /= INSULATED) THEN
      E_WALLB = 0._EB
      VOLSUM  = 0._EB
      DO N=1,ONE_D%N_MATL
         IF (ONE_D%MATL_COMP(N)%RHO(NWP)<=TWO_EPSILON_EB) CYCLE
         ML => MATERIAL(ONE_D%MATL_INDEX(N))
         VOLSUM  = VOLSUM  + ONE_D%MATL_COMP(N)%RHO(NWP)/RHO_ADJUSTED(LAYER_INDEX(NWP),N)
         E_WALLB = E_WALLB + ONE_D%MATL_COMP(N)%RHO(NWP)*ML%EMISSIVITY/RHO_ADJUSTED(LAYER_INDEX(NWP),N)
      ENDDO
      IF (VOLSUM > 0._EB) E_WALLB = E_WALLB/VOLSUM
   ENDIF

   ! Get heat losses from convection and radiation out of back of surface

   SELECT CASE(BACKING)

      CASE(VOID)  ! Non-insulated backing to an ambient void

         IF (SF%TMP_GAS_BACK>0._EB) THEN
            TMP_GAS_BACK = TMP_0(BC%KK) + EVALUATE_RAMP(T-T_BEGIN,SF%RAMP(TIME_TGB)%INDEX)*(SF%TMP_GAS_BACK-TMP_0(BC%KK))
         ELSE
            TMP_GAS_BACK = TMP_0(BC%KK)
         ENDIF
         DTMP = TMP_GAS_BACK - B1%TMP_B
         HTCB = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC_B,SF)
         Q_CON_B = HTCB*DTMP
         Q_RAD_IN_B   =  E_WALLB*SIGMA*TMP_GAS_BACK**4
         Q_LIQUID_B = 0._EB
         LAYER_DIVIDE = REAL(SF%N_LAYERS+1)
         MF_FRAC = 1._EB

      CASE(INSULATED)  ! No heat transfer out the back

         HTCB      = 0._EB
         Q_CON_B   = 0._EB
         Q_RAD_IN_B   = 0._EB
         E_WALLB   = 0._EB
         Q_LIQUID_B = 0._EB
         TMP_GAS_BACK = B1%TMP_B

      CASE(EXPOSED)  ! The backside is exposed to gas in current or adjacent mesh.

         Q_LIQUID_B = 0._EB

         IF (BACK_MESH/=NM .AND. BACK_MESH>0) THEN  ! Back side is in other mesh.
            TMP_GAS_BACK = B1_BACK%TMP_G
            DTMP = TMP_GAS_BACK - B1%TMP_B
            IF (.NOT.ISOLATED_THIN_WALL_BACK) THEN
               HTCB = HEAT_TRANSFER_COEFFICIENT(BACK_MESH,DTMP,SF_HTC_B,SF_BACK,WALL_INDEX_IN=BACK_WALL_INDEX)
            ELSE
               HTCB = 0._EB
            ENDIF
            Q_RAD_IN_B  = B1_BACK%Q_RAD_IN
         ELSE  ! Back side is in current mesh.
            TMP_GAS_BACK  = TMP(BC_BACK%IIG,BC_BACK%JJG,BC_BACK%KKG)
            DTMP = TMP_GAS_BACK - B1%TMP_B
            IF (PRESENT(WALL_INDEX)) &
               HTCB = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC_B,SF_BACK,WALL_INDEX_IN=ONE_D%BACK_INDEX)
            IF (PRESENT(THIN_WALL_INDEX)) THEN
               IF (BACK_WALL_INDEX>0) THEN
                  HTCB = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC_B,SF_BACK,WALL_INDEX_IN=BACK_WALL_INDEX)
               ELSE
                  HTCB = 0._EB
               ENDIF
            ENDIF
            IF (PRESENT(CFACE_INDEX)) &
               HTCB = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC_B,SF_BACK,CFACE_INDEX_IN=ONE_D%BACK_INDEX)
            B1_BACK%HEAT_TRANS_COEF = HTCB
            Q_RAD_IN_B  = B1_BACK%Q_RAD_IN
            IF (N_LP_ARRAY_INDICES>0 .AND. PRESENT(WALL_INDEX)) Q_LIQUID_B = -SUM(B2_BACK%LP_CPUA(:)) + B1%Q_CONDENSE
         ENDIF
         Q_CON_B = HTCB*DTMP

   END SELECT

   ! Get total thickness of solid and compute radius for cylindrical and spherical coordinate systems.

   THICKNESS = SUM(ONE_D%LAYER_THICKNESS(1:ONE_D%N_LAYERS))

   IF (SF%GEOMETRY==SURF_INNER_CYLINDRICAL) THEN
      DO I=0,NWP
         R_S(I) = SF%INNER_RADIUS + ONE_D%X(I)
      ENDDO
   ELSE
      DO I=0,NWP
         R_S(I) = SF%INNER_RADIUS + ONE_D%X(NWP) - ONE_D%X(I)
      ENDDO
   ENDIF

   ! Calculate reaction rates based on the solid phase reactions

   Q_S = 0._EB

   PYROLYSIS_PREDICTED_IF: IF (ONE_D%PYROLYSIS_MODEL==PYROLYSIS_PREDICTED) THEN

      T_BOIL_EFF = TMPA
      CALL PERFORM_PYROLYSIS

   ELSEIF (ONE_D%PYROLYSIS_MODEL==PYROLYSIS_SPECIFIED) THEN PYROLYSIS_PREDICTED_IF

      ! Take off energy corresponding to specified burning rate
      !***** Figure out what to do for adjust if different spec have different HOC

      Q_S(1) = Q_S(1) - B1%M_DOT_G_PP_ADJUST(REACTION(1)%FUEL_SMIX_INDEX)*SF%H_V/DX_S(1)

   ENDIF PYROLYSIS_PREDICTED_IF

   ! Determine additional heat sources
   
   Q_ADD = 0._EB
   Q_IR = 0._EB

   ! Add internal heat source specified by user

   DO I=1,NWP
      Q_ADD(I) = ONE_D%HEAT_SOURCE(LAYER_INDEX(I))*EVALUATE_RAMP(T-T_BEGIN,ONE_D%RAMP_IHS_INDEX(LAYER_INDEX(I)))
   ENDDO

   ! Add special convection term for Boundary Fuel Model

   IF (SF%BOUNDARY_FUEL_MODEL) THEN
      B1%HEAT_TRANS_COEF = 0._EB
      HTCB = 0._EB
      Q_CON_F = 0._EB
      Q_CON_B = 0._EB
      N = 0
      DO I=1,NWP
         IF (SF%SURFACE_VOLUME_RATIO(LAYER_INDEX(I))<=0._EB) CYCLE
         DTMP = B1%TMP_G - ONE_D%TMP(I)
         IF (PRESENT(WALL_INDEX)) THEN
            HTCF = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC_F,SF,WALL_INDEX_IN=WALL_INDEX)
         ELSEIF(PRESENT(CFACE_INDEX)) THEN
            HTCF = HEAT_TRANSFER_COEFFICIENT(NM,DTMP,SF_HTC_F,SF,CFACE_INDEX_IN=CFACE_INDEX)
         ENDIF
         DEL_DOT_Q_SC = HTCF*DTMP
         Q_ADD(I) = Q_ADD(I) + SF%SURFACE_VOLUME_RATIO(LAYER_INDEX(I))*SF%PACKING_RATIO(LAYER_INDEX(I))*DEL_DOT_Q_SC
         ! Track average h_c for computing h_m in SURFACE_OXIDATION_MODEL
         B1%HEAT_TRANS_COEF = B1%HEAT_TRANS_COEF + HTCF
         N = N + 1
      ENDDO
      B1%HEAT_TRANS_COEF = B1%HEAT_TRANS_COEF/REAL(N,EB)
      HTCF = 0._EB
   ENDIF

   ! Calculate internal radiation for Cartesian geometry only

   IF (ONE_D%INTERNAL_RADIATION) THEN
      DO I=1,NWP
         IF (SF%KAPPA_S(LAYER_INDEX(I))<0._EB) THEN
            VOLSUM = 0._EB
            KAPSUM = 0._EB
            DO N=1,ONE_D%N_MATL
               IF (ONE_D%MATL_COMP(N)%RHO(I)<=TWO_EPSILON_EB) CYCLE
               ML  => MATERIAL(ONE_D%MATL_INDEX(N))
               VOLSUM = VOLSUM + ONE_D%MATL_COMP(N)%RHO(I)/RHO_ADJUSTED(LAYER_INDEX(I),N)
               KAPSUM = KAPSUM + ONE_D%MATL_COMP(N)%RHO(I)*ML%KAPPA_S/RHO_ADJUSTED(LAYER_INDEX(I),N)
            ENDDO
            IF (VOLSUM>0._EB) TWO_DX_KAPPA_S(I) = 2._EB*KAPSUM/(RDX_S(I)*VOLSUM)
         ELSE
            TWO_DX_KAPPA_S(I) = 2._EB*SF%KAPPA_S(LAYER_INDEX(I))/RDX_S(I)
         ENDIF
      ENDDO
      ! solution inwards
      RFLUX_UP = B1%Q_RAD_IN + (1._EB-B1%EMISSIVITY)*Q_RAD_OUT_OLD/(B1%EMISSIVITY+1.0E-10_EB)
      DO I=1,NWP
         RFLUX_DOWN =  ( RFLUX_UP + TWO_DX_KAPPA_S(I)*SIGMA*ONE_D%TMP(I)**4 ) / (1._EB + TWO_DX_KAPPA_S(I))
         Q_IR(I) = Q_IR(I) + (RFLUX_UP - RFLUX_DOWN)*RDX_S(I)
         RFLUX_UP = RFLUX_DOWN
      ENDDO
      ! solution outwards
      RFLUX_UP = Q_RAD_IN_B + (1._EB-E_WALLB)*RFLUX_UP
      DO I=NWP,1,-1
         RFLUX_DOWN =  ( RFLUX_UP + TWO_DX_KAPPA_S(I)*SIGMA*ONE_D%TMP(I)**4 ) / (1._EB + TWO_DX_KAPPA_S(I))
         Q_IR(I) = Q_IR(I) + (RFLUX_UP - RFLUX_DOWN)*RDX_S(I)
         RFLUX_UP = RFLUX_DOWN
      ENDDO
      Q_RAD_OUT = B1%EMISSIVITY*RFLUX_DOWN
   ENDIF

   ! Add internal radiation and additional heat sources to pyrolysis
   Q_S = Q_S + Q_IR + Q_ADD
   ! If the 3D solver is used, divide Q_S by 3
   Q_S = Q_S/REAL(SF%HT_DIM,EB)

   ! Explicitly update the temperature field and adjust time step if the change in temperature exceeds DELTA_TMP_MAX

   IF (ICYC>WALL_INCREMENT) THEN
      IF (ONE_D%INTERNAL_RADIATION) THEN
         Q_NET_F = Q_CON_F
         Q_NET_B = Q_CON_B
      ELSE
         Q_NET_F = B1%Q_RAD_IN - B1%EMISSIVITY*SIGMA*B1%TMP_F**4 + Q_CON_F
         Q_NET_B = Q_RAD_IN_B  - E_WALLB      *SIGMA*B1%TMP_B**4 + Q_CON_B
      ENDIF
      DO I=1,NWP
         DELTA_TMP(I) = DT_BC*Q_S(I)/ONE_D%RHO_C_S(I)
      ENDDO
      TMP_RATIO = MAX(TWO_EPSILON_EB,MAXVAL(ABS(DELTA_TMP(1:NWP)))/SF%DELTA_TMP_MAX)
      DT_BC_SUB_OLD = DT_BC_SUB
      DT_BC_SUB = DT_BC/REAL(MIN(NINT(SF%TIME_STEP_FACTOR*WALL_INCREMENT),MAX(1,NINT(TMP_RATIO))),EB)
      DT_BC_SUB = MIN( DT_BC-T_BC_SUB , DT_BC_SUB , DT_FO )
      ! If DT change, rebuild Q_S
      IF (ONE_D%PYROLYSIS_MODEL==PYROLYSIS_PREDICTED .AND. DT_BC_SUB_OLD/=DT_BC_SUB) THEN
         Q_S = 0._EB
         CALL PERFORM_PYROLYSIS
         ! Add internal radiation and additional heat sources to pyrolysis
         Q_S = Q_S + Q_IR + Q_ADD
         ! If the 3D solver is used, divide Q_S by 3
         Q_S = Q_S/REAL(SF%HT_DIM,EB)
      ENDIF
   ENDIF

   T_BC_SUB = T_BC_SUB + DT_BC_SUB

   ! Store the mass and energy fluxes from this time sub-iteration

   IF (ONE_D%INTERNAL_RADIATION) THEN
      Q_RAD_OUT_OLD = Q_RAD_OUT
      B1%Q_RAD_OUT = B1%Q_RAD_OUT + Q_RAD_OUT*DT_BC_SUB/DT_BC
   ENDIF

   IF (ONE_D%PYROLYSIS_MODEL==PYROLYSIS_PREDICTED .AND. .NOT.PRESENT(THIN_WALL_INDEX)) THEN
      RDT = 1._EB/(DT_BC*REAL(SF%HT_DIM,EB))
      IF (PRESENT(WALL_INDEX)) THEN
         IF (OBSTRUCTION(OBST_INDEX)%THIN .AND. .NOT.CELL(CELL_INDEX(BC%II,BC%JJ,BC%KK))%SOLID) RDT = 1._EB/DT_BC
      ENDIF
      Q_DOT_G_PP_NET                 = Q_DOT_G_PP_NET             + Q_DOT_G_PP*DT_BC_SUB*RDT
      Q_DOT_O2_PP_NET                = Q_DOT_O2_PP_NET            + Q_DOT_O2_PP*DT_BC_SUB*RDT
      M_DOT_G_PP_ADJUST_NET          = M_DOT_G_PP_ADJUST_NET      + B1%AREA_ADJUST*M_DOT_G_PP_ADJUST*DT_BC_SUB*RDT
      M_DOT_G_PP_ACTUAL_NET          = M_DOT_G_PP_ACTUAL_NET      +                M_DOT_G_PP_ACTUAL*DT_BC_SUB*RDT
      M_DOT_S_PP_NET(1:ONE_D%N_MATL) = M_DOT_S_PP_NET(1:ONE_D%N_MATL) + M_DOT_S_PP(1:ONE_D%N_MATL)*DT_BC_SUB*RDT
   ENDIF

   ! Adjust the material layer masses and thicknesses

   REMESH_LAYER = .FALSE.
   TMP_CHECK = .FALSE.

   PYROLYSIS_PREDICTED_IF_2: IF (ONE_D%PYROLYSIS_MODEL==PYROLYSIS_PREDICTED) THEN

      ! Convert Q_S to kW
      DO I=1,NWP
         Q_S(I) = Q_S(I)*ABS(R_S(I-1)**I_GRAD-R_S(I)**I_GRAD)
      ENDDO

      POINT_LOOP2: DO I=1,NWP
         
         NODE_RDT(I) = 1.E9_EB

         ! No need for a computation if no reactions
         
         IF (ALL(ABS(RHO_DOT(:,I))<TWO_EPSILON_EB)) THEN
            REGRID_FACTOR(I) = 1._EB   
            CYCLE POINT_LOOP2 
         ELSE
            REGRID_FACTOR(I) = 0._EB
         ENDIF

         ! Compute regrid factors

         MATERIAL_LOOP1a: DO N=1,ONE_D%N_MATL
            ONE_D%MATL_COMP(N)%RHO(I) = MAX( 0._EB , ONE_D%MATL_COMP(N)%RHO(I) - DT_BC_SUB*RHO_DOT(N,I) )
            IF (RHO_DOT(N,I) > TWO_EPSILON_EB) NODE_RDT(I) = MIN(NODE_RDT(I),MATERIAL(ONE_D%MATL_INDEX(N))%RENODE_DELTA_T)
            REGRID_FACTOR(I) = REGRID_FACTOR(I) + ONE_D%MATL_COMP(N)%RHO(I)/RHO_ADJUSTED(LAYER_INDEX(I),N)
         ENDDO MATERIAL_LOOP1a

         ! If there is any non-shrinking material, the material matrix will remain, and no shrinking is allowed

         MATERIAL_LOOP1b: DO N=1,ONE_D%N_MATL
            IF (ONE_D%MATL_COMP(N)%RHO(I)<=TWO_EPSILON_EB) CYCLE MATERIAL_LOOP1b
            ML  => MATERIAL(ONE_D%MATL_INDEX(N))
            IF (.NOT. ML%ALLOW_SHRINKING) THEN
               REGRID_FACTOR(I) = MAX(REGRID_FACTOR(I),1._EB)
               EXIT MATERIAL_LOOP1b
            ENDIF
         ENDDO MATERIAL_LOOP1b

         ! If there is any non-swelling material, the material matrix will remain, and no swelling is allowed

         MATERIAL_LOOP1c: DO N=1,ONE_D%N_MATL
            IF (ONE_D%MATL_COMP(N)%RHO(I)<=TWO_EPSILON_EB) CYCLE MATERIAL_LOOP1c
            ML  => MATERIAL(ONE_D%MATL_INDEX(N))
            IF (.NOT. ML%ALLOW_SWELLING) THEN
               REGRID_FACTOR(I) = MIN(REGRID_FACTOR(I),1._EB)
               EXIT MATERIAL_LOOP1c
            ENDIF
         ENDDO MATERIAL_LOOP1c

         ! In points that change thickness, update the density

         IF (ABS(REGRID_FACTOR(I)-1._EB)>=TWO_EPSILON_EB) THEN
            MATERIAL_LOOP1d: DO N=1,ONE_D%N_MATL
               IF(REGRID_FACTOR(I)>TWO_EPSILON_EB) ONE_D%MATL_COMP(N)%RHO(I) = ONE_D%MATL_COMP(N)%RHO(I)/REGRID_FACTOR(I)
            ENDDO MATERIAL_LOOP1d
         ENDIF

      ENDDO POINT_LOOP2

      ! Compute new coordinates if the solid changes thickness. Save new coordinates in X_S_NEW.
      ! Remesh layer if any node goes to zero thickness

      R_S_NEW(NWP) = 0._EB
      DO I=NWP-1,0,-1
         R_S_NEW(I) = ( R_S_NEW(I+1)**I_GRAD + (R_S(I)**I_GRAD-R_S(I+1)**I_GRAD)*REGRID_FACTOR(I+1) )**(1./REAL(I_GRAD,EB))
      ENDDO

      X_S_NEW(0) = 0._EB
      CELL_ZERO = .FALSE.

      DO I=1,NWP
         X_S_NEW(I) = R_S_NEW(0) - R_S_NEW(I)
         ! If Cell disappears we must remesh
         IF ((X_S_NEW(I)-X_S_NEW(I-1)) < TWO_EPSILON_EB) THEN
            REMESH_LAYER(LAYER_INDEX(I)) = .TRUE.
            CELL_ZERO(LAYER_INDEX(I)) = .TRUE.
         ELSE
            ! If cell size changes enough compared to prior remseh size, remesh
            IF (ABS((X_S_NEW(I)-X_S_NEW(I-1))/ONE_D%DX_OLD(I)-1._EB) > SF%REMESH_RATIO) REMESH_LAYER(LAYER_INDEX(I)) = .TRUE.
            IF (I > 1 .AND. REMESH_LAYER(LAYER_INDEX(I))) THEN
               IF (LAYER_INDEX(I-1)==LAYER_INDEX(I) .AND. ABS(ONE_D%TMP(I)-ONE_D%TMP(I-1)) > NODE_RDT(I)) THEN
                  IF (.NOT. CELL_ZERO(LAYER_INDEX(I))) REMESH_LAYER(LAYER_INDEX(I)) = .FALSE.
                  TMP_CHECK(LAYER_INDEX(I)) = .TRUE.
               ENDIF
            ENDIF
            IF (I < NWP .AND. REMESH_LAYER(LAYER_INDEX(I))) THEN
               IF (LAYER_INDEX(I+1)==LAYER_INDEX(I) .AND. ABS(ONE_D%TMP(I+1)-ONE_D%TMP(I)) > NODE_RDT(I)) THEN
                  IF (.NOT. CELL_ZERO(LAYER_INDEX(I))) REMESH_LAYER(LAYER_INDEX(I)) = .FALSE.
                  TMP_CHECK(LAYER_INDEX(I)) = .TRUE.
               ENDIF
            ENDIF               
         ENDIF
      ENDDO

      ! Check for layers that are too small, layer thickness dropping too much, or couldn't drop enough nodes last remesh
      I = 0
      DO NL=1,ONE_D%N_LAYERS
         ONE_D%LAYER_THICKNESS(NL) = X_S_NEW(I+ONE_D%N_LAYER_CELLS(NL)) - X_S_NEW(I)
         IF (ONE_D%N_LAYER_CELLS(NL) > 0 .AND. ONE_D%LAYER_THICKNESS(NL) < 0.1_EB*ONE_D%MINIMUM_LAYER_THICKNESS(NL)) THEN
            REMESH_LAYER(NL) = .TRUE.
         ELSE
            IF (.NOT. TMP_CHECK(NL) .AND. ONE_D%LAYER_THICKNESS_OLD(NL)-ONE_D%LAYER_THICKNESS(NL) > &
               1.5_EB*MINVAL(ONE_D%DX_OLD(I+1:I+ONE_D%N_LAYER_CELLS(NL)))) THEN
                  REMESH_LAYER(NL) = .TRUE.
            ELSE
               IF (.NOT. TMP_CHECK(NL) .AND. ONE_D%REMESH_NWP(NL) > ONE_D%N_LAYER_CELLS(NL)) REMESH_LAYER(NL) = .TRUE.
            ENDIF
         ENDIF
         I = I + ONE_D%N_LAYER_CELLS(NL)
      ENDDO
      
      ! If any nodes go to zero, apportion Q_S to surrounding nodes.

      IF (ANY(CELL_ZERO(1:ONE_D%N_LAYERS)) .AND. NWP > 1) THEN
         IF (X_S_NEW(1)-X_S_NEW(0) < TWO_EPSILON_EB) Q_S(2) = Q_S(2) + Q_S(1)
         IF (X_S_NEW(NWP)-X_S_NEW(NWP-1) < TWO_EPSILON_EB) Q_S(NWP-1) = Q_S(NWP-1) + Q_S(NWP)
         DO I=2,NWP-1
            IF (X_S_NEW(I) - X_S_NEW(I-1) < TWO_EPSILON_EB) THEN
               N = 0
               IF (X_S_NEW(I-1) - X_S_NEW(I-2) > TWO_EPSILON_EB) N=N+1
               IF (X_S_NEW(I+1) - X_S_NEW(I) > TWO_EPSILON_EB) N=N+2
               SELECT CASE (N)
                  CASE(1)
                     Q_S(I-1) = Q_S(I-1) + Q_S(I)
                  CASE(2)
                     Q_S(I+1) = Q_S(I+1) + Q_S(I)
                  CASE(3)
                     VOL = (R_S_NEW(I-1)**I_GRAD-R_S_NEW(I)**I_GRAD) / &
                           ((R_S_NEW(I-1)**I_GRAD-R_S_NEW(I)**I_GRAD)+(R_S_NEW(I)**I_GRAD-R_S_NEW(I+1)**I_GRAD))
                     Q_S(I-1) = Q_S(I-1) + Q_S(I) * VOL
                     Q_S(I+1) = Q_S(I+1) + Q_S(I) * (1._EB-VOL)
               END SELECT
            ENDIF
         ENDDO
      ENDIF

      REMESH_CHECK = ANY(ABS(REGRID_FACTOR(1:NWP)-1._EB)>TWO_EPSILON_EB)

      ! Some node changes size but no layer trips remesh check. Just redo node weight with X_S_NEW      
     
      IF (REMESH_CHECK .AND. ALL(.NOT. REMESH_LAYER(1:ONE_D%N_LAYERS))) THEN
         CALL GET_WALL_NODE_WEIGHTS(NWP,ONE_D%N_LAYERS,ONE_D%N_LAYER_CELLS(1:ONE_D%N_LAYERS), &
                                   ONE_D%LAYER_THICKNESS(1:ONE_D%N_LAYERS),SF%GEOMETRY,X_S_NEW(0:NWP),LAYER_DIVIDE,DX_S(1:NWP), &
                                   RDX_S(0:NWP+1),RDXN_S(0:NWP),DX_WGT_S(0:NWP),DXF,DXB,LAYER_INDEX(0:NWP+1),MF_FRAC(1:NWP),&
                                   SF%INNER_RADIUS)
      ENDIF

      N_LAYER_CELLS_NEW(1:ONE_D%N_LAYERS) = ONE_D%N_LAYER_CELLS(1:ONE_D%N_LAYERS)      

      ! Some layer needs to be checked for a remesh.

      REMESH_LAYER_1: IF (ANY(REMESH_LAYER(1:ONE_D%N_LAYERS))) THEN
         NWP_NEW = 0
         THICKNESS = 0._EB
         I = 0
         DX_MIN = 0._EB
         ONE_D%X(0:NWP) = X_S_NEW(0:NWP)
      
         LAYER_LOOP: DO NL=1,ONE_D%N_LAYERS

            ! Layer is too small. Delete it and shift any following nodes up in x-distance

            IF (ONE_D%LAYER_THICKNESS(NL) < 0.1_EB*ONE_D%MINIMUM_LAYER_THICKNESS(NL)) THEN
               N_LAYER_CELLS_NEW(NL) = 0
               ONE_D%X(I+ONE_D%N_LAYER_CELLS(NL):NWP) = ONE_D%X(I+ONE_D%N_LAYER_CELLS(NL):NWP)-ONE_D%LAYER_THICKNESS(NL)
               ONE_D%DX_OLD(I+1:I+ONE_D%N_LAYER_CELLS(NL)) = 0._EB
               ONE_D%LAYER_THICKNESS(NL) = 0._EB
               ONE_D%LAYER_THICKNESS_OLD(NL) = 0._EB
               ONE_D%REMESH_NWP(NL) = 0
               CYCLE LAYER_LOOP
            ENDIF

            ! Layer is one cell. Only invoke remeshing if the layer grows in size which could add a wall node
            
            IF (ONE_D%N_LAYER_CELLS(NL) == 1) THEN
               ! Don't remesh layer if one cell and shrinking
               IF (REMESH_LAYER(NL) .AND. REGRID_FACTOR(I+1) < 1._EB) REMESH_LAYER(NL) = .FALSE.
               IF (.NOT. REMESH_LAYER(NL)) THEN
                  N_LAYER_CELLS_NEW(NL) = 1
                  NWP_NEW = NWP_NEW + 1
                  I = I + 1
                  THICKNESS = THICKNESS + ONE_D%LAYER_THICKNESS(NL)
                  ONE_D%SMALLEST_CELL_SIZE(NL) = ONE_D%LAYER_THICKNESS(NL)
                  CYCLE LAYER_LOOP
               ENDIF
            ENDIF

            ! Layer isn't tagged for possible remeshing. Just update new cell count and thickness.

            IF (.NOT. REMESH_LAYER(NL)) THEN
               N_LAYER_CELLS_NEW(NL) = ONE_D%N_LAYER_CELLS(NL)
               NWP_NEW = NWP_NEW + N_LAYER_CELLS_NEW(NL)
               DO N = I,I+ONE_D%N_LAYER_CELLS(NL)-1
                  IF (ONE_D%X(I+1)-ONE_D%X(I) < ONE_D%SMALLEST_CELL_SIZE(NL)) ONE_D%SMALLEST_CELL_SIZE(NL) = ONE_D%X(I+1)-ONE_D%X(I)
               ENDDO
               THICKNESS = THICKNESS + ONE_D%LAYER_THICKNESS(NL)
               I = I + ONE_D%N_LAYER_CELLS(NL)
               CYCLE LAYER_LOOP
            ENDIF
      
            ! Layer passes all checks for a possible remesh

            CALL GET_N_LAYER_CELLS(ONE_D%MIN_DIFFUSIVITY(NL),ONE_D%LAYER_THICKNESS(NL),ONE_D%STRETCH_FACTOR(NL),&
                                   ONE_D%CELL_SIZE_FACTOR(NL),ONE_D%CELL_SIZE(NL),ONE_D%N_LAYER_CELLS_MAX(NL),&
                                   N_LAYER_CELLS_NEW(NL),ONE_D%SMALLEST_CELL_SIZE(NL),ONE_D%DDSUM(NL))
            ONE_D%REMESH_NWP(NL) = N_LAYER_CELLS_NEW(NL)
            LAYER_CELL_CHECK: IF (ONE_D%N_LAYER_CELLS(NL) - N_LAYER_CELLS_NEW(NL) > 1) THEN
               N_LAYER_CELLS_NEW(NL) = ONE_D%N_LAYER_CELLS(NL)- 1
               IF (MOD(N_LAYER_CELLS_NEW(NL),2)==0) THEN
                  DDSUM = 0._EB
                  DO N=1,N_LAYER_CELLS_NEW(NL)/2
                     DDSUM = DDSUM + ONE_D%STRETCH_FACTOR(NL)**(N-1)
                  ENDDO
                  DDSUM = 2._EB*DDSUM
               ELSE
                  DDSUM = 0._EB
                  DO N=1,(N_LAYER_CELLS_NEW(NL)-1)/2
                     DDSUM = DDSUM + ONE_D%STRETCH_FACTOR(NL)**(N-1)
                  ENDDO
                  DDSUM = 2._EB*DDSUM + ONE_D%STRETCH_FACTOR(NL)**((N_LAYER_CELLS_NEW(NL)-1)/2)
               ENDIF
               ONE_D%SMALLEST_CELL_SIZE(NL) = ONE_D%LAYER_THICKNESS(NL) / DDSUM
               ONE_D%DDSUM(NL) = DDSUM
            ENDIF LAYER_CELL_CHECK
            
            THICKNESS = THICKNESS + ONE_D%LAYER_THICKNESS(NL)
            I = I + N_LAYER_CELLS_NEW(NL)       
            NWP_NEW = NWP_NEW + N_LAYER_CELLS_NEW(NL)

         ENDDO LAYER_LOOP

         IF (NWP_NEW > ONE_D%N_CELLS_MAX) THEN
            WRITE(MESSAGE,'(A,I5,A,A)') 'ERROR(300): N_LAYER_CELLS_MAX should be at least ',NWP_NEW,' for ',TRIM(SF%ID)
            CALL SHUTDOWN(MESSAGE,PROCESS_0_ONLY=.FALSE.)
         ENDIF
                
      ELSE

         THICKNESS = SUM(ONE_D%LAYER_THICKNESS)
         ONE_D%X(0:NWP) = X_S_NEW(0:NWP)
         
      ENDIF REMESH_LAYER_1

      ! Shrinking wall has gone to zero thickness.

      IF (THICKNESS<=TWO_EPSILON_EB) THEN
         ONE_D%TMP(0:NWP+1)                    = MAX(TMPMIN,TMP_GAS_BACK)
         B1%TMP_F                              = MIN(TMPMAX,MAX(TMPMIN,TMP_GAS_BACK))
         B1%TMP_B                              = MIN(TMPMAX,MAX(TMPMIN,TMP_GAS_BACK))
         B1%Q_CON_F                            = 0._EB
         B1%Q_RAD_OUT                          = 0._EB
         ONE_D%N_LAYER_CELLS(1:ONE_D%N_LAYERS) = 0
         B1%BURNAWAY       = .TRUE.
         ONE_D%PART_MASS(1:ONE_D%N_LPC) = 0._EB
         ONE_D%PART_ENTHALPY(1:ONE_D%N_LPC) = 0._EB
         B1%T_MATL_PART = 0._EB
         B1%M_DOT_PART_ACTUAL = 0._EB
         IF (OBST_INDEX > 0) THEN
            IF (OBSTRUCTION(OBST_INDEX)%CONSUMABLE) OBSTRUCTION(OBST_INDEX)%MASS = -1.
         ENDIF
         EXIT SUB_TIMESTEP_LOOP
      ENDIF      
      
      ! Re-generate grid for a wall changing thickness
      
      REMESH_LAYER_2: IF (ANY(REMESH_LAYER(1:ONE_D%N_LAYERS))) THEN
         RHO_H_S = 0._EB
         TMP_S = 0._EB

         ! Store wall enthalpy for later temperature extraction.

         DO I=1,NWP
            VOL = (THICKNESS+SF%INNER_RADIUS-ONE_D%X(I-1))**I_GRAD-(THICKNESS+SF%INNER_RADIUS-ONE_D%X(I))**I_GRAD
            MATL_REMESH: DO N=1,ONE_D%N_MATL
               IF (ONE_D%MATL_COMP(N)%RHO(I)<=TWO_EPSILON_EB) CYCLE MATL_REMESH
               ML  => MATERIAL(ONE_D%MATL_INDEX(N))
               ITMP = MIN(I_MAX_TEMP-1,INT(ONE_D%TMP(I)))
               H_S = ML%H(ITMP)+(ONE_D%TMP(I)-REAL(ITMP,EB))*(ML%H(ITMP+1)-ML%H(ITMP))
               RHO_H_S(I) = RHO_H_S(I) + ONE_D%MATL_COMP(N)%RHO(I) * H_S
            ENDDO MATL_REMESH
            RHO_H_S(I) = RHO_H_S(I) * VOL
            DO N=1,ONE_D%N_MATL
               ONE_D%MATL_COMP(N)%RHO(I) = ONE_D%MATL_COMP(N)%RHO(I) * VOL
            ENDDO
            TMP_S(I)=ONE_D%TMP(I)*VOL
         ENDDO

         CALL GET_WALL_NODE_COORDINATES(NWP_NEW,NWP,ONE_D%N_LAYERS,N_LAYER_CELLS_NEW(1:ONE_D%N_LAYERS), &
            ONE_D%N_LAYER_CELLS(1:ONE_D%N_LAYERS),ONE_D%SMALLEST_CELL_SIZE(1:ONE_D%N_LAYERS), &
            ONE_D%STRETCH_FACTOR(1:ONE_D%N_LAYERS),REMESH_LAYER(1:ONE_D%N_LAYERS),X_S_NEW(0:NWP_NEW),ONE_D%X(0:NWP), &
            ONE_D%LAYER_THICKNESS(1:ONE_D%N_LAYERS))
         CALL GET_WALL_NODE_WEIGHTS(NWP_NEW,ONE_D%N_LAYERS,N_LAYER_CELLS_NEW(1:ONE_D%N_LAYERS),ONE_D%LAYER_THICKNESS,SF%GEOMETRY, &
            X_S_NEW(0:NWP_NEW),LAYER_DIVIDE,DX_S(1:NWP_NEW),RDX_S(0:NWP_NEW+1),RDXN_S(0:NWP_NEW),&
            DX_WGT_S(0:NWP_NEW),DXF,DXB,LAYER_INDEX(0:NWP_NEW+1),MF_FRAC(1:NWP_NEW),SF%INNER_RADIUS)
         IF(NWP_NEW < NWP) ONE_D%DX_OLD(NWP_NEW:NWP) = 0._EB ! Zero out old values if needed 
         ONE_D%DX_OLD(1:NWP_NEW) = DX_S(1:NWP_NEW)
         ONE_D%LAYER_THICKNESS_OLD(1:ONE_D%N_LAYERS) = ONE_D%LAYER_THICKNESS(1:ONE_D%N_LAYERS)

         ! Interpolate densities and temperature from old grid to new grid

         ALLOCATE(INT_WGT(NWP_NEW,NWP),STAT=IZERO)
         CALL GET_INTERPOLATION_WEIGHTS(SF%GEOMETRY,NWP,NWP_NEW,SF%INNER_RADIUS,ONE_D%X(0:NWP),X_S_NEW(0:NWP_NEW),INT_WGT)
         N_CELLS = MAX(NWP,NWP_NEW)
         CALL INTERPOLATE_WALL_ARRAY(N_CELLS,NWP,NWP_NEW,INT_WGT,Q_S(1:N_CELLS))
         CALL INTERPOLATE_WALL_ARRAY(N_CELLS,NWP,NWP_NEW,INT_WGT,RHO_H_S(1:N_CELLS))
         CALL INTERPOLATE_WALL_ARRAY(N_CELLS,NWP,NWP_NEW,INT_WGT,TMP_S(1:N_CELLS))

         DO I=1,NWP_NEW
            VOL = (THICKNESS+SF%INNER_RADIUS-X_S_NEW(I-1))**I_GRAD-(THICKNESS+SF%INNER_RADIUS-X_S_NEW(I))**I_GRAD
            TMP_S(I) = TMP_S(I) / VOL
         ENDDO

         DO N=1,ONE_D%N_MATL
            ML  => MATERIAL(ONE_D%MATL_INDEX(N))
            CALL INTERPOLATE_WALL_ARRAY(N_CELLS,NWP,NWP_NEW,INT_WGT,ONE_D%MATL_COMP(N)%RHO(1:N_CELLS))
         ENDDO

         DEALLOCATE(INT_WGT)

         ! Extract temperature

         DO I=1,NWP_NEW
            H_NODE = RHO_H_S(I)
            T_NODE = TMP_S(I)
            ITER = 0
            T_SEARCH: DO
               ITER = ITER + 1
               C_S = 0._EB
               H_S = 0._EB
               CP1 = 0
               CP2 = 0
               ITMP = MIN(I_MAX_TEMP-1,INT(T_NODE))
               H_S = 0._EB
               T_S: DO N=1,ONE_D%N_MATL
                  IF (ONE_D%MATL_COMP(N)%RHO(I)<=0._EB) CYCLE T_S
                  ML  => MATERIAL(ONE_D%MATL_INDEX(N))
                  H_S = H_S + (ML%H(ITMP)+(T_NODE-REAL(ITMP,EB))*(ML%H(ITMP+1)-ML%H(ITMP)))*ONE_D%MATL_COMP(N)%RHO(I)
                  CP1 = CP1 + ML%H(ITMP)/REAL(ITMP,EB)*ONE_D%MATL_COMP(N)%RHO(I)
                  CP2 = CP2 + ML%H(ITMP+1)/REAL(ITMP+1,EB)*ONE_D%MATL_COMP(N)%RHO(I)
               ENDDO T_S
               C_S = H_S/T_NODE
               DENOM = C_S+T_NODE*(CP2-CP1)
               IF (ABS(DENOM) < TWO_EPSILON_EB) THEN
                  ONE_D%TMP(I) = T_NODE
               ELSE
                  ONE_D%TMP(I) = T_NODE + (H_NODE - H_S)/DENOM
               ENDIF
               IF (ABS(ONE_D%TMP(I) - T_NODE) < 0.0001_EB) EXIT T_SEARCH
               IF (ITER > 20) THEN
                  ONE_D%TMP(I) = 0.5_EB*(ONE_D%TMP(I)+T_NODE)
                  EXIT T_SEARCH
               ENDIF
               T_NODE = ONE_D%TMP(I)
            ENDDO T_SEARCH
            DO N=1,ONE_D%N_MATL
               ONE_D%MATL_COMP(N)%RHO(I) = ONE_D%MATL_COMP(N)%RHO(I) /&
                  ((SF%INNER_RADIUS+X_S_NEW(NWP_NEW)-X_S_NEW(I-1))**I_GRAD-(SF%INNER_RADIUS+X_S_NEW(NWP_NEW)-X_S_NEW(I))**I_GRAD)
            ENDDO
         ENDDO

         ONE_D%TMP(0)         = 2._EB*B1%TMP_F-ONE_D%TMP(1)       ! Make sure front surface temperature stays the same
         ONE_D%TMP(NWP_NEW+1) = 2._EB*B1%TMP_B-ONE_D%TMP(NWP_NEW) ! Make sure back surface temperature stays the same

         ONE_D%N_LAYER_CELLS(1:ONE_D%N_LAYERS) = N_LAYER_CELLS_NEW(1:ONE_D%N_LAYERS)
         NWP = NWP_NEW
         ONE_D%X(0:NWP) = X_S_NEW(0:NWP)      ! Note: X(NWP+1...) are not set to zero.
      ENDIF REMESH_LAYER_2

      ! Convert Q_S back to kW/m^3

      DO I=1,NWP
         Q_S(I) = Q_S(I)/((SF%INNER_RADIUS+ONE_D%X(NWP)-ONE_D%X(I-1))**I_GRAD-(SF%INNER_RADIUS+ONE_D%X(NWP)-ONE_D%X(I))**I_GRAD)
      ENDDO

   ENDIF PYROLYSIS_PREDICTED_IF_2

   ! Calculate thermal properties

   ONE_D%K_S = 0._EB
   RHO_S   = 0._EB
   ONE_D%RHO_C_S = 0._EB
   B1%EMISSIVITY = 0._EB
   E_FOUND = .FALSE.

   POINT_LOOP3: DO I=1,NWP
      VOLSUM = 0._EB
      ITMP = MIN(I_MAX_TEMP-1,INT(ONE_D%TMP(I)))
      MATERIAL_LOOP3: DO N=1,ONE_D%N_MATL
         IF (ONE_D%MATL_COMP(N)%RHO(I)<=TWO_EPSILON_EB) CYCLE MATERIAL_LOOP3
         ML  => MATERIAL(ONE_D%MATL_INDEX(N))
         VOLSUM = VOLSUM + ONE_D%MATL_COMP(N)%RHO(I)/RHO_ADJUSTED(LAYER_INDEX(I),N)
         ONE_D%K_S(I) = ONE_D%K_S(I) + ONE_D%MATL_COMP(N)%RHO(I)*ML%K_S(ITMP)/RHO_ADJUSTED(LAYER_INDEX(I),N)
         ONE_D%RHO_C_S(I) = ONE_D%RHO_C_S(I) + ONE_D%MATL_COMP(N)%RHO(I)*ML%C_S(ITMP)

         IF (.NOT.E_FOUND) B1%EMISSIVITY = B1%EMISSIVITY + ONE_D%MATL_COMP(N)%RHO(I)*ML%EMISSIVITY/RHO_ADJUSTED(LAYER_INDEX(I),N)
         RHO_S(I) = RHO_S(I) + ONE_D%MATL_COMP(N)%RHO(I)
      ENDDO MATERIAL_LOOP3
      IF (SF%PACKING_RATIO(LAYER_INDEX(I))>0._EB) ONE_D%K_S(I) = ONE_D%K_S(I)*SF%PACKING_RATIO(LAYER_INDEX(I))

      IF (VOLSUM > 0._EB) THEN
         ONE_D%K_S(I) = ONE_D%K_S(I)/VOLSUM
         IF (.NOT.E_FOUND) B1%EMISSIVITY = B1%EMISSIVITY/VOLSUM
      ENDIF
      IF (B1%EMISSIVITY>=0._EB) E_FOUND = .TRUE.
      IF (ONE_D%K_S(I)<=TWO_EPSILON_EB)      ONE_D%K_S(I)      = 10000._EB
      IF (ONE_D%RHO_C_S(I)<=TWO_EPSILON_EB)  ONE_D%RHO_C_S(I)  = 0.001_EB
   ENDDO POINT_LOOP3

   IF (SF%EMISSIVITY_SPECIFIED) B1%EMISSIVITY = SF%EMISSIVITY

   ! Calculate average K_S at grid cell boundaries. Store result in K_S

   ONE_D%K_S(0)     = ONE_D%K_S(1)
   ONE_D%K_S(NWP+1) = ONE_D%K_S(NWP)
   DO I=1,NWP-1
      ONE_D%K_S(I)  = 1._EB / ( DX_WGT_S(I)/ONE_D%K_S(I) + (1._EB-DX_WGT_S(I))/ONE_D%K_S(I+1) )
   ENDDO

   ! Update the 1-D heat transfer equation

   KODXF = ONE_D%K_S(0)/DXF
   KODXB = ONE_D%K_S(NWP)/DXB

   DO I=1,NWP
      BBS(I) = -0.5_EB*DT_BC_SUB*ONE_D%K_S(I-1)*RDXN_S(I-1)*RDX_S(I)/ONE_D%RHO_C_S(I)
      AAS(I) = -0.5_EB*DT_BC_SUB*ONE_D%K_S(I)  *RDXN_S(I)  *RDX_S(I)/ONE_D%RHO_C_S(I)
   ENDDO
   DDS(1:NWP) = 1._EB - AAS(1:NWP) - BBS(1:NWP)
   DO I=1,NWP
      CCS(I) = ONE_D%TMP(I) - AAS(I)*(ONE_D%TMP(I+1)-ONE_D%TMP(I)) + BBS(I)*(ONE_D%TMP(I)-ONE_D%TMP(I-1)) &
               + DT_BC_SUB*Q_S(I)/ONE_D%RHO_C_S(I)
   ENDDO

   IF (SF%DIRICHLET_FRONT) THEN
      RFACF2 = -1._EB
      QDXKF  = 2._EB*B1%TMP_F
   ELSEIF (.NOT.RADIATION .OR. ONE_D%INTERNAL_RADIATION .OR. ISOLATED_THIN_WALL) THEN
      RFACF  = 0.5_EB*HTCF
      RFACF2 = (KODXF-RFACF)/(KODXF+RFACF)
      QDXKF  = (HTCF*B1%TMP_G  + Q_LIQUID_F)/(KODXF+RFACF)
   ELSE
      RFACF  = 0.5_EB*HTCF + 2._EB*B1%EMISSIVITY*SIGMA*B1%TMP_F**3
      RFACF2 = (KODXF-RFACF)/(KODXF+RFACF)
      QDXKF  = (HTCF*B1%TMP_G  + Q_LIQUID_F + B1%Q_RAD_IN + 3._EB*B1%EMISSIVITY*SIGMA*B1%TMP_F**4) / (KODXF+RFACF)
   ENDIF

   IF (SF%DIRICHLET_BACK) THEN
      RFACB2 = -1._EB
      QDXKB  = 2._EB*B1%TMP_B
   ELSEIF (.NOT.RADIATION .OR. ONE_D%INTERNAL_RADIATION .OR. ISOLATED_THIN_WALL_BACK) THEN
      RFACB  = 0.5_EB*HTCB
      RFACB2 = (KODXB-RFACB)/(KODXB+RFACB)
      QDXKB  = (HTCB*TMP_GAS_BACK + Q_LIQUID_B)/(KODXB+RFACB)
   ELSE
      RFACB  = 0.5_EB*HTCB + 2._EB*E_WALLB*SIGMA*B1%TMP_B**3
      RFACB2 = (KODXB-RFACB)/(KODXB+RFACB)
      QDXKB  = (HTCB*TMP_GAS_BACK + Q_LIQUID_B + Q_RAD_IN_B + 3._EB*E_WALLB*SIGMA*B1%TMP_B**4) / (KODXB+RFACB)
   ENDIF

   CCS(1)   = CCS(1)   - BBS(1)  *QDXKF
   CCS(NWP) = CCS(NWP) - AAS(NWP)*QDXKB
   DDT(1:NWP) = DDS(1:NWP)
   DDT(1)   = DDT(1)   + BBS(1)  *RFACF2
   DDT(NWP) = DDT(NWP) + AAS(NWP)*RFACB2

   TRIDIAGONAL_SOLVER_1: DO I=2,NWP
      RR     = BBS(I)/DDT(I-1)
      DDT(I) = DDT(I) - RR*AAS(I-1)
      CCS(I) = CCS(I) - RR*CCS(I-1)
   ENDDO TRIDIAGONAL_SOLVER_1
   CCS(NWP)  = CCS(NWP)/DDT(NWP)
   TRIDIAGONAL_SOLVER_2: DO I=NWP-1,1,-1
      CCS(I) = (CCS(I) - AAS(I)*CCS(I+1))/DDT(I)
   ENDDO TRIDIAGONAL_SOLVER_2

   DELTA_TMP(1:NWP) = MIN(TMPMAX,MAX(TMPMIN,CCS(1:NWP))) - ONE_D%TMP(1:NWP)
   ONE_D%TMP(1:NWP) = ONE_D%TMP(1:NWP) + DELTA_TMP(1:NWP)
   DELTA_TMP(0)     = MAX(TMPMIN,ONE_D%TMP(1)  *RFACF2+QDXKF) - ONE_D%TMP(0) ! Ghost values, allow them to be large
   ONE_D%TMP(0)     = ONE_D%TMP(0)     + DELTA_TMP(0)
   DELTA_TMP(NWP+1) = MAX(TMPMIN,ONE_D%TMP(NWP)*RFACB2+QDXKB) - ONE_D%TMP(NWP+1)
   ONE_D%TMP(NWP+1) = ONE_D%TMP(NWP+1) + DELTA_TMP(NWP+1)
   ONE_D%DELTA_TMP(0:NWP+1) = ONE_D%DELTA_TMP(0:NWP+1) + DELTA_TMP(0:NWP+1)

   B1%Q_CON_F = B1%Q_CON_F + HTCF*DT_BC_SUB*(B1%TMP_G-0.5_EB*B1%TMP_F)
   B1%TMP_F_OLD = B1%TMP_F  ! Save this value for output of effective HTC

   B1%TMP_F  = 0.5_EB*(ONE_D%TMP(0)+ONE_D%TMP(1))
   B1%TMP_B  = 0.5_EB*(ONE_D%TMP(NWP)+ONE_D%TMP(NWP+1))

   B1%Q_CON_F = B1%Q_CON_F - 0.5_EB*HTCF*DT_BC_SUB*B1%TMP_F
   IF (RADIATION .AND. .NOT.ONE_D%INTERNAL_RADIATION) B1%Q_RAD_OUT = B1%Q_RAD_OUT + &
                                                      DT_BC_SUB*(B1%TMP_F_OLD**4+2._EB*B1%TMP_F_OLD**3*(B1%TMP_F-B1%TMP_F_OLD))

   ! Clipping for excessively high or low temperatures

   B1%TMP_F  = MIN(TMPMAX,MAX(TMPMIN,B1%TMP_F))
   B1%TMP_B  = MIN(TMPMAX,MAX(TMPMIN,B1%TMP_B))

   ! Updated particle production

   IF (ONE_D%N_LPC > 0) THEN
      ONE_D%PART_MASS(1:ONE_D%N_LPC) = ONE_D%PART_MASS(1:ONE_D%N_LPC) + DT_BC_SUB * M_DOT_PART_S(1:ONE_D%N_LPC)
      ONE_D%PART_ENTHALPY(1:ONE_D%N_LPC) = ONE_D%PART_ENTHALPY(1:ONE_D%N_LPC) + DT_BC_SUB * Q_DOT_PART_S(1:ONE_D%N_LPC)
      B1%T_MATL_PART = B1%T_MATL_PART + DT_BC_SUB
      B1%M_DOT_PART_ACTUAL = SUM(M_DOT_PART_S(1:ONE_D%N_LPC))
   ENDIF

   ! Determine if the iterations are done, otherwise return to the top

   IF (T_BC_SUB>=DT_BC-TWO_EPSILON_EB) EXIT SUB_TIMESTEP_LOOP

   B1%N_SUBSTEPS = B1%N_SUBSTEPS + 1

ENDDO SUB_TIMESTEP_LOOP

IF (ALLOCATED(B1%M_DOT_G_PP_ACTUAL) .AND. ONE_D%PYROLYSIS_MODEL==PYROLYSIS_PREDICTED) THEN
   B1%M_DOT_G_PP_ADJUST             = M_DOT_G_PP_ADJUST_NET
   B1%M_DOT_G_PP_ACTUAL             = M_DOT_G_PP_ACTUAL_NET
   B1%Q_DOT_G_PP                    = Q_DOT_G_PP_NET
   B1%Q_DOT_O2_PP                   = Q_DOT_O2_PP_NET
   ONE_D%M_DOT_S_PP(1:ONE_D%N_MATL) = M_DOT_S_PP_NET(1:ONE_D%N_MATL)
ENDIF
IF (ALLOCATED(RHO_DOT) .AND. ONE_D%PYROLYSIS_MODEL==PYROLYSIS_PREDICTED) DEALLOCATE(RHO_DOT)

B1%Q_CON_F = B1%Q_CON_F / DT_BC
IF (RADIATION .AND. .NOT. ONE_D%INTERNAL_RADIATION) B1%Q_RAD_OUT = B1%Q_RAD_OUT / DT_BC * SIGMA * B1%EMISSIVITY
IF (.NOT. SF%BOUNDARY_FUEL_MODEL) B1%HEAT_TRANS_COEF = HTCF

! If any gas massflux or particle mass flux is non-zero or the surface temperature exceeds the ignition temperature,
! set the ignition time

IF (B1%T_IGN > T) THEN
   IF (SUM(B1%M_DOT_G_PP_ADJUST(1:N_TRACKED_SPECIES)) > 0._EB .OR. B1%M_DOT_PART_ACTUAL > 0._EB ) B1%T_IGN = T
   IF (B1%TMP_F>=SF%TMP_IGN) B1%T_IGN = T
ENDIF

! If the surface temperature is less than the extinction temperature, stop the burning

IF (SF%TMP_IGN<50000._EB .AND. B1%TMP_F<SF%TMP_EXT .AND. B1%T_IGN<T) B1%T_IGN = HUGE(1._EB)

CONTAINS


!> \brief Compute the solid phase reaction term for all solid phase nodes

SUBROUTINE PERFORM_PYROLYSIS

USE PHYSICAL_FUNCTIONS, ONLY: GET_MASS_FRACTION,GET_SPECIFIC_HEAT
REAL(EB), DIMENSION(N_TRACKED_SPECIES) :: M_DOT_G_PPP_ADJUST,M_DOT_G_PPP_ACTUAL,ZZ_GET
REAL(EB), DIMENSION(MAX_MATERIALS) :: RHO_TEMP,M_DOT_S_PPP
REAL(EB) :: Q_DOT_G_PPP,Q_DOT_O2_PPP,CP_FILM,TMP_FILM,H_MASS,&
            Y_O2_G,Y_O2_F,M_DOT_O2_PP,Y_LOWER,Y_UPPER,M_DOT_ERROR,M_DOT_ERROR_OLD,Y_O2_F_OLD,DY,DE
REAL(EB), DIMENSION(MAX_LPC) :: Q_DOT_PART,M_DOT_PART
INTEGER :: ITER,MAX_ITER
LOGICAL :: REMOVE_LAYER
REAL(EB), PARAMETER :: M_DOT_ERROR_TOL=1.E-6_EB

! Get surface oxygen mass fraction

IF (O2_INDEX>0) THEN
   ZZ_GET(1:N_TRACKED_SPECIES) = MAX(0._EB,ZZ(BC%IIG,BC%JJG,BC%KKG,1:N_TRACKED_SPECIES))
   CALL GET_MASS_FRACTION(ZZ_GET,O2_INDEX,Y_O2_F)
ELSE
   Y_O2_F = 0._EB
ENDIF

MAX_ITER=1

! Initialize parameters for oxidative pyrolysis mass transfer model

IF (OXPYRO_MODEL) THEN
   MAX_ITER=20
   Y_O2_G  = Y_O2_F
   Y_LOWER = 0._EB
   Y_UPPER = Y_O2_G
   Y_O2_F  = MIN(B2%Y_O2_F,Y_O2_G) ! initial guess for Newton method (B2 only available for wall cell right now)
   TMP_FILM = (B1%TMP_F+TMP(BC%IIG,BC%JJG,BC%KKG))/2._EB
   CALL GET_SPECIFIC_HEAT(ZZ_GET,CP_FILM,TMP_FILM)
   H_MASS = B1%HEAT_TRANS_COEF/CP_FILM
   M_DOT_ERROR = 0._EB
   DY = 0._EB
ENDIF

O2_LOOP: DO ITER=1,MAX_ITER

   ! Set mass and energy fluxes to zero for this time sub-iteration

   M_DOT_G_PP_ADJUST(1:N_TRACKED_SPECIES) = 0._EB
   M_DOT_G_PP_ACTUAL(1:N_TRACKED_SPECIES) = 0._EB
   M_DOT_S_PP(1:ONE_D%N_MATL)             = 0._EB
   Q_DOT_G_PP                             = 0._EB
   Q_DOT_O2_PP                            = 0._EB
   M_DOT_PART_S                           = 0._EB
   Q_DOT_PART_S                           = 0._EB

   ! Loop over all solid cells and compute the reaction rate of each material component, N, in each cell, I, RHO_DOT(N,I)

   POINT_LOOP1: DO I=1,NWP

      ! Create a temporary array to hold the material component densities at the current depth layer, I

      DO N=1,ONE_D%N_MATL
         RHO_TEMP(N) = ONE_D%MATL_COMP(N)%RHO(I)
      ENDDO

      IF (ONE_D%LAYER_THICKNESS(LAYER_INDEX(I))<ONE_D%MINIMUM_LAYER_THICKNESS(LAYER_INDEX(I))) THEN
         REMOVE_LAYER = .TRUE.
         B1%LAYER_REMOVED = .TRUE.
      ELSE
         REMOVE_LAYER = .FALSE.
      ENDIF

      IF (PRESENT(PARTICLE_INDEX)) THEN
         CALL PYROLYSIS(ONE_D%N_MATL,ONE_D%MATL_INDEX,SURF_INDEX,BC%IIG,BC%JJG,BC%KKG,ONE_D%TMP(I),B1%TMP_F,Y_O2_F,BC%IOR,&
                        RHO_DOT(1:ONE_D%N_MATL,I),RHO_TEMP(1:ONE_D%N_MATL),ONE_D%X(I-1),DX_S,DT_BC_SUB,&
                        M_DOT_G_PPP_ADJUST,M_DOT_G_PPP_ACTUAL,M_DOT_S_PPP,Q_S(I),Q_DOT_G_PPP,Q_DOT_O2_PPP,&
                        Q_DOT_PART,M_DOT_PART,T_BOIL_EFF,B1%B_NUMBER,LAYER_INDEX(I),REMOVE_LAYER,ONE_D,B1,SOLID_CELL_INDEX=I,&
                        R_DROP=R_SURF,LPU=U_SURF,LPV=V_SURF,LPW=W_SURF)
      ELSE
         CALL PYROLYSIS(ONE_D%N_MATL,ONE_D%MATL_INDEX,SURF_INDEX,BC%IIG,BC%JJG,BC%KKG,ONE_D%TMP(I),B1%TMP_F,Y_O2_F,BC%IOR,&
                        RHO_DOT(1:ONE_D%N_MATL,I),RHO_TEMP(1:ONE_D%N_MATL),ONE_D%X(I-1),DX_S,DT_BC_SUB,&
                        M_DOT_G_PPP_ADJUST,M_DOT_G_PPP_ACTUAL,M_DOT_S_PPP,Q_S(I),Q_DOT_G_PPP,Q_DOT_O2_PPP,&
                        Q_DOT_PART,M_DOT_PART,T_BOIL_EFF,B1%B_NUMBER,LAYER_INDEX(I),REMOVE_LAYER,ONE_D,B1,SOLID_CELL_INDEX=I)
      ENDIF

      ! Sum the mass and heat generation within the solid layers (PPP) and transfer result to the surface (PP).
      ! GEOM_FACTOR accounts for cylindrical and spherical geometry. If this is not a Lagrangian particle, scale with
      ! the original radius (SF%THICKNESS) because the VENT to which the flux is applied does not change in area.

      IF (ONE_D%N_LAYERS==1 .AND. REMOVE_LAYER) MF_FRAC(I) = 1._EB

      IF (PRESENT(PARTICLE_INDEX)) THEN
         GEOM_FACTOR = MF_FRAC(I)*(R_S(I-1)**I_GRAD-R_S(I)**I_GRAD)/(I_GRAD*R_S(0)**(I_GRAD-1))
      ELSE
         GEOM_FACTOR = MF_FRAC(I)*(R_S(I-1)**I_GRAD-R_S(I)**I_GRAD)/(I_GRAD*(SF%THICKNESS+SF%INNER_RADIUS)**(I_GRAD-1))
      ENDIF

      Q_DOT_G_PP  = Q_DOT_G_PP  + Q_DOT_G_PPP *GEOM_FACTOR
      Q_DOT_O2_PP = Q_DOT_O2_PP + Q_DOT_O2_PPP*GEOM_FACTOR

      M_DOT_G_PP_ADJUST = M_DOT_G_PP_ADJUST + M_DOT_G_PPP_ADJUST*GEOM_FACTOR
      M_DOT_G_PP_ACTUAL = M_DOT_G_PP_ACTUAL + M_DOT_G_PPP_ACTUAL*GEOM_FACTOR

      M_DOT_S_PP(1:ONE_D%N_MATL) = M_DOT_S_PP(1:ONE_D%N_MATL)  + M_DOT_S_PPP(1:ONE_D%N_MATL)*GEOM_FACTOR

      ! Compute particle mass flux at the surface

      IF (ONE_D%N_LPC > 0) THEN
         M_DOT_PART_S(1:ONE_D%N_LPC) = M_DOT_PART_S(1:ONE_D%N_LPC) + GEOM_FACTOR * M_DOT_PART(1:ONE_D%N_LPC)
         Q_DOT_PART_S(1:ONE_D%N_LPC) = Q_DOT_PART_S(1:ONE_D%N_LPC) + GEOM_FACTOR * Q_DOT_PART(1:ONE_D%N_LPC)
      ENDIF

   ENDDO POINT_LOOP1

   !======== Mass transfer resistance to surface O2 concentration =============
   OXPYRO_MODEL_IF: IF (OXPYRO_MODEL) THEN
      ! Solve for Y_O2_F using bounded Newton method
      M_DOT_ERROR_OLD = M_DOT_ERROR
      M_DOT_O2_PP = DOT_PRODUCT(Z2Y(O2_INDEX,1:N_TRACKED_SPECIES),M_DOT_G_PP_ACTUAL(1:N_TRACKED_SPECIES))
      M_DOT_ERROR = M_DOT_O2_PP - H_MASS*(Y_O2_F - Y_O2_G)

      IF (ABS(M_DOT_ERROR)<M_DOT_ERROR_TOL) THEN
         B2%Y_O2_F = Y_O2_F
         B2%Y_O2_ITER = ITER
         EXIT O2_LOOP
      ENDIF

      IF (M_DOT_ERROR<0._EB) THEN
         Y_UPPER = Y_O2_F
      ELSE
         Y_LOWER = Y_O2_F
      ENDIF
      Y_O2_F_OLD = Y_O2_F
      DE = M_DOT_ERROR-M_DOT_ERROR_OLD
      IF (ABS(DE)>TWO_EPSILON_EB) THEN
         ! use Newton
         Y_O2_F = MIN(Y_UPPER,MAX(Y_LOWER, Y_O2_F_OLD - M_DOT_ERROR*DY/DE ))
      ELSE
         ! bisect
         Y_O2_F = 0.5_EB*(Y_LOWER+Y_UPPER)
      ENDIF
      DY = Y_O2_F - Y_O2_F_OLD
   ENDIF OXPYRO_MODEL_IF
   !===========================================================================

ENDDO O2_LOOP

END SUBROUTINE PERFORM_PYROLYSIS

END SUBROUTINE SOLID_HEAT_TRANSFER


!> \brief Calculate the solid phase reaction. Return heat and mass generation rates per unit volume.
!>
!> \param N_MATS Number of material components in the solid
!> \param MATL_INDEX (1:N_MATS) Indices of the material components from the master material list
!> \param SURF_INDEX Index of surface, used only for liquids
!> \param IIG I index of nearest gas phase cell
!> \param JJG J index of nearest gas phase cell
!> \param KKG K index of nearest gas phase cell
!> \param TMP_S Solid interior temperature (K)
!> \param TMP_F Solid surface temperature (K)
!> \param Y_O2_F Solid surface value of oxygen mass fraction
!> \param IOR Index of orientation of the surface with the liquid droplet, if appropropriate (0 for gas phase droplet)
!> \param RHO_DOT_OUT (1:N_MATS) Array of component reaction rates (kg/m3/s)
!> \param RHO_S (1:N_MATS) Array of component densities (kg/m3)
!> \param DEPTH Distance from surface (m)
!> \param DX_S Array of node sizes (m)
!> \param DT_BC Time step used by the solid phase solver (s)
!> \param M_DOT_G_PPP_ADJUST (1:N_TRACKED_SPECIES) Adjusted mass generation rate per unit volume of the gas species
!> \param M_DOT_G_PPP_ACTUAL (1:N_TRACKED_SPECIES) Actual mass generation rate per unit volume of the gas species
!> \param M_DOT_S_PPP (1:N_MATS) Mass generation/depletion rate per unit volume of solid components (kg/m3/s)
!> \param Q_DOT_S_PPP Heat release rate per unit volume (W/m3)
!> \param Q_DOT_G_PPP Rate of energy required to bring gaseous pyrolyzate to the surrounding gas temperature (W/m3)
!> \param Q_DOT_O2_PPP Heat release rate per unit volume due to char oxidation in grid cell abutting surface (W/m3)
!> \param Q_DOT_PART Rate of enthalpy production of particles created in reactions (J/m3/s)
!> \param M_DOT_PART Rate of mass production of particles created in reactions (kg/m3/s)
!> \param T_BOIL_EFF Effective boiling temperature (K)
!> \param B_NUMBER B-number of liquid surface
!> \param LAYER_INDEX Index of the solid layer
!> \param REMOVE_LAYER Indicates if the entire layer is being removed
!> \param ONE_D Pointer to collection of variables for the interior of the solid
!> \param B1 Pointer to surface variables
!> \param SOLID_CELL_INDEX (OPTIONAL) Index of the interior solid cell
!> \param R_DROP (OPTIONAL) Radius of liquid droplet
!> \param LPU (OPTIONAL) x component of droplet velocity (m/s)
!> \param LPV (OPTIONAL) y component of droplet velocity (m/s)
!> \param LPW (OPTIONAL) z component of droplet velocity (m/s)

SUBROUTINE PYROLYSIS(N_MATS,MATL_INDEX,SURF_INDEX,IIG,JJG,KKG,TMP_S,TMP_F,Y_O2_F,IOR,RHO_DOT_OUT,RHO_S,DEPTH,DX_S,DT_BC,&
                     M_DOT_G_PPP_ADJUST,M_DOT_G_PPP_ACTUAL,M_DOT_S_PPP,Q_DOT_S_PPP,Q_DOT_G_PPP,Q_DOT_O2_PPP,&
                     Q_DOT_PART,M_DOT_PART,T_BOIL_EFF,B_NUMBER,LAYER_INDEX,REMOVE_LAYER,ONE_D,B1,SOLID_CELL_INDEX,&
                     R_DROP,LPU,LPV,LPW)

USE PHYSICAL_FUNCTIONS, ONLY: GET_MASS_FRACTION,GET_VISCOSITY,GET_PARTICLE_ENTHALPY,GET_SPECIFIC_HEAT,&
                              GET_MASS_FRACTION_ALL,GET_EQUIL_DATA,GET_SENSIBLE_ENTHALPY,GET_Y_SURF,GET_FILM_PROPERTIES
USE MATH_FUNCTIONS, ONLY: INTERPOLATE1D_UNIFORM
USE TURBULENCE, ONLY: RAYLEIGH_HEAT_FLUX_MODEL,RAYLEIGH_MASS_FLUX_MODEL
INTEGER, INTENT(IN) :: N_MATS,SURF_INDEX,IIG,JJG,KKG,IOR,LAYER_INDEX
INTEGER, INTENT(IN), OPTIONAL :: SOLID_CELL_INDEX
LOGICAL, INTENT(IN) :: REMOVE_LAYER
REAL(EB), INTENT(OUT), DIMENSION(:,:) :: RHO_DOT_OUT(N_MATS)
REAL(EB), INTENT(IN) :: TMP_S,TMP_F,DT_BC,DEPTH,RHO_S(N_MATS),Y_O2_F
REAL(EB), INTENT(IN), OPTIONAL :: R_DROP,LPU,LPV,LPW
REAL(EB), INTENT(IN), DIMENSION(NWP_MAX) :: DX_S
REAL(EB), DIMENSION(:) :: ZZ_GET(1:N_TRACKED_SPECIES),Y_ALL(1:N_SPECIES)
REAL(EB), DIMENSION(:), INTENT(OUT) :: M_DOT_G_PPP_ADJUST(N_TRACKED_SPECIES),M_DOT_G_PPP_ACTUAL(N_TRACKED_SPECIES)
REAL(EB), DIMENSION(:), INTENT(OUT) :: M_DOT_S_PPP(MAX_MATERIALS),Q_DOT_PART(MAX_LPC),M_DOT_PART(MAX_LPC)
REAL(EB), INTENT(OUT) :: Q_DOT_S_PPP,Q_DOT_G_PPP,Q_DOT_O2_PPP,B_NUMBER
REAL(EB), INTENT(INOUT) :: T_BOIL_EFF(MAX_MATERIALS)
INTEGER, INTENT(IN), DIMENSION(:) :: MATL_INDEX(N_MATS)
INTEGER :: N,NN,NNN,J,NS,SMIX_INDEX(N_MATS),NWP,NP,NP2,ITMP
TYPE(MATERIAL_TYPE), POINTER :: ML
TYPE(SURFACE_TYPE), POINTER :: SF
TYPE(BOUNDARY_ONE_D_TYPE), POINTER :: ONE_D
TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1
REAL(EB) :: REACTION_RATE,Y_O2,X_O2,MW(N_MATS),Y_GAS(N_MATS),Y_TMP(N_MATS),Y_SV(N_MATS),X_SV(N_MATS),X_L(N_MATS),&
            D_FILM,H_MASS,RE_L,SHERWOOD,MFLUX,MU_FILM,SC_FILM,TMP_FILM,TMP_G,U2,V2,W2,VEL,&
            DR,R_S_0,R_S_1,H_R,H_R_B,H_S_B,H_S,LENGTH_SCALE,SUM_Y_GAS,SUM_Y_SV,NU_O2_CHAR,Y_O2_S,&
            SUM_Y_SV_SMIX(N_TRACKED_SPECIES),X_L_SUM,RHO_DOT_EXTRA,MFLUX_MAX,RHO_FILM,CP_FILM,PR_FILM,K_FILM,EVAP_FILM_FAC,&
            RHO_DOT,RHO_DOT_REAC(MAX_REACTIONS),RHO_DOT_REAC_SUM,H_MASS_DNS
LOGICAL :: LIQUID(N_MATS),SPEC_ID_ALREADY_USED(N_MATS),DO_EVAPORATION

B_NUMBER = 0._EB
Q_DOT_S_PPP = 0._EB
Q_DOT_G_PPP = 0._EB
Q_DOT_O2_PPP = 0._EB
M_DOT_S_PPP = 0._EB
M_DOT_G_PPP_ADJUST = 0._EB
M_DOT_G_PPP_ACTUAL = 0._EB
M_DOT_PART = 0._EB
Q_DOT_PART = 0._EB
RHO_DOT_OUT = 0._EB
SF => SURFACE(SURF_INDEX)

! Determine if any liquids are present. If they are, determine if this is a the surface layer.

DO_EVAPORATION = .FALSE.
IF (ANY(MATERIAL(MATL_INDEX(:))%PYROLYSIS_MODEL==PYROLYSIS_LIQUID)) THEN
   IF (PRESENT(SOLID_CELL_INDEX)) THEN
      IF (SOLID_CELL_INDEX==1) DO_EVAPORATION = .TRUE.
   ENDIF
ENDIF

! If this is surface liquid layer, calculate the Spalding B number and other liquid-specific variables

IF_DO_EVAPORATION: IF (DO_EVAPORATION) THEN

   ! Calculate a sum needed to calculate the volume fraction of liquid components

   LIQUID  = .FALSE.
   X_L_SUM = 0._EB
   MATERIAL_LOOP_00: DO N=1,N_MATS
      ML => MATERIAL(MATL_INDEX(N))
      IF (ML%PYROLYSIS_MODEL/=PYROLYSIS_LIQUID) CYCLE MATERIAL_LOOP_00
      IF (RHO_S(N) < TWO_EPSILON_EB) CYCLE MATERIAL_LOOP_00
      LIQUID(N) = .TRUE.
      X_L_SUM = X_L_SUM + RHO_S(N)/ML%RHO_S
   ENDDO MATERIAL_LOOP_00

   IF (X_L_SUM < TWO_EPSILON_EB) RETURN

   SUM_Y_GAS = 0._EB
   SPEC_ID_ALREADY_USED = .FALSE.
   SMIX_INDEX = 0
   X_L = 0._EB
   X_SV = 0._EB

   MATERIAL_LOOP_0: DO N=1,N_MATS

      IF (.NOT.LIQUID(N)) CYCLE MATERIAL_LOOP_0
      ML => MATERIAL(MATL_INDEX(N))

      SMIX_INDEX(N) = MAXLOC(ML%NU_GAS(:,1),1)
      ZZ_GET(1:N_TRACKED_SPECIES) = MAX(0._EB,ZZ(IIG,JJG,KKG,1:N_TRACKED_SPECIES))

      IF (ML%MW<0._EB) THEN  ! No molecular weight specified; assume the liquid component evaporates into the defined gas species
         MW(N) = SPECIES_MIXTURE(SMIX_INDEX(N))%MW
      ELSE  ! the user has specified a molecular weight for the liquid component
         MW(N) = ML%MW
      ENDIF

      ! Determine the mass fraction of evaporated MATL N in the first gas phase grid cell

      IF (SPECIES_MIXTURE(SMIX_INDEX(N))%SINGLE_SPEC_INDEX > 0) THEN
         CALL GET_MASS_FRACTION_ALL(ZZ_GET,Y_ALL)
         Y_GAS(N) = Y_ALL(SPECIES_MIXTURE(SMIX_INDEX(N))%SINGLE_SPEC_INDEX)
         IF (SPECIES_MIXTURE(SMIX_INDEX(N))%CONDENSATION_SMIX_INDEX > 0) &
               Y_GAS(N) = Y_GAS(N) - ZZ_GET(SPECIES_MIXTURE(SMIX_INDEX(N))%CONDENSATION_SMIX_INDEX)
      ELSE
         Y_GAS(N) = ZZ_GET(SMIX_INDEX(N))
      ENDIF

      ! Determine volume fraction of MATL N in the liquid and then the surface vapor layer

      T_BOIL_EFF(N) = ML%TMP_BOIL
      CALL GET_EQUIL_DATA(MW(N),TMP_F,PBAR(KKG,PRESSURE_ZONE(IIG,JJG,KKG)),H_R,H_R_B,T_BOIL_EFF(N),X_SV(N),ML%H_R(1,:))
      X_L(N)  = RHO_S(N)/(ML%RHO_S*X_L_SUM)  ! Volume fraction of MATL component N in the liquid
      X_SV(N) = X_L(N)*X_SV(N)               ! Volume fraction of MATL component N in the surface vapor based on Raoult's law

      ! Calculate sums to be used to compute B number
      IF (.NOT.SPEC_ID_ALREADY_USED(N)) SUM_Y_GAS = SUM_Y_GAS + Y_GAS(N)
      SPEC_ID_ALREADY_USED(N) = .TRUE.

   ENDDO MATERIAL_LOOP_0

   ! Convert mole fraction to mass fraction
   CALL GET_Y_SURF(N_MATS,ZZ_GET,X_SV,Y_SV,MW,SMIX_INDEX)

   ! Compute the Spalding B number

   SUM_Y_SV = SUM(Y_SV)
   SUM_Y_SV_SMIX = 0._EB
   MATERIAL_LOOP_1: DO N=1,N_MATS
      IF (.NOT.LIQUID(N)) CYCLE MATERIAL_LOOP_1
      SUM_Y_SV_SMIX(SMIX_INDEX(N)) = SUM_Y_SV_SMIX(SMIX_INDEX(N)) + Y_SV(N)
   ENDDO MATERIAL_LOOP_1

   IF (SUM_Y_SV<ONE_M_EPS) THEN
      B_NUMBER = MAX(0._EB,(SUM_Y_SV-SUM(Y_GAS))/(1._EB-SUM_Y_SV))
   ELSE
      B_NUMBER = 1.E6_EB  ! Fictitiously high B number intended to push mass flux to its upper limit
   ENDIF

   ! Compute an effective gas phase mass fraction, Y_GAS, corresponding to each liquid component, N

   Y_TMP = 0._EB
   MATERIAL_LOOP_2: DO N=1,N_MATS
      IF (.NOT.LIQUID(N)) CYCLE MATERIAL_LOOP_2
      IF (SUM_Y_SV_SMIX(SMIX_INDEX(N))>TWO_EPSILON_EB) Y_TMP(N) = Y_SV(N)*Y_GAS(N)/SUM_Y_SV_SMIX(SMIX_INDEX(N))
   ENDDO MATERIAL_LOOP_2
   Y_GAS = Y_TMP

   ! Get film properties

   SELECT CASE(SF%GEOMETRY)
      CASE DEFAULT                                  ; EVAP_FILM_FAC = PLATE_FILM_FAC
      CASE(SURF_SPHERICAL)                          ; EVAP_FILM_FAC = SPHERE_FILM_FAC
      CASE(SURF_CYLINDRICAL,SURF_INNER_CYLINDRICAL) ; EVAP_FILM_FAC = PLATE_FILM_FAC
      CASE(SURF_CARTESIAN)                          ; EVAP_FILM_FAC = PLATE_FILM_FAC
   END SELECT

   CALL GET_FILM_PROPERTIES(N_MATS,EVAP_FILM_FAC,Y_SV,Y_GAS,SMIX_INDEX,TMP_F,TMP(IIG,JJG,KKG),ZZ_GET,&
                           PBAR(KKG,PRESSURE_ZONE(IIG,JJG,KKG)),TMP_FILM,MU_FILM,K_FILM,CP_FILM,D_FILM,&
                           RHO_FILM,PR_FILM,SC_FILM)

   ! Compute mass transfer coefficient

   H_MASS_IF: IF (SF%HM_FIXED>=0._EB) THEN

      H_MASS = SF%HM_FIXED

   ELSE H_MASS_IF

      SELECT CASE(ABS(IOR))
         CASE(0); H_MASS_DNS = 0._EB
         CASE(1); H_MASS_DNS = 2._EB*D_FILM*RDX(IIG)
         CASE(2); H_MASS_DNS = 2._EB*D_FILM*RDY(JJG)
         CASE(3); H_MASS_DNS = 2._EB*D_FILM*RDZ(KKG)
      END SELECT
      
      IF (SIM_MODE==DNS_MODE) THEN 

         H_MASS = H_MASS_DNS

      ELSE

         IF (PRESENT(LPU) .AND. PRESENT(LPV) .AND. PRESENT(LPW)) THEN
            U2 = 0.5_EB*(U(IIG,JJG,KKG)+U(IIG-1,JJG,KKG))
            V2 = 0.5_EB*(V(IIG,JJG,KKG)+V(IIG,JJG-1,KKG))
            W2 = 0.5_EB*(W(IIG,JJG,KKG)+W(IIG,JJG,KKG-1))
            VEL = SQRT((U2-LPU)**2+(V2-LPV)**2+(W2-LPW)**2)
         ELSE
            VEL = SQRT(2._EB*KRES(IIG,JJG,KKG))
         ENDIF
         CALL GET_VISCOSITY(ZZ_GET,MU_FILM,TMP_FILM)
         IF (PRESENT(R_DROP)) THEN
            LENGTH_SCALE = 2._EB*R_DROP
         ELSE

            LENGTH_SCALE = SF%CONV_LENGTH
         ENDIF
         RE_L     = RHO_FILM*VEL*LENGTH_SCALE/MU_FILM
         SELECT CASE(SF%GEOMETRY)
            CASE DEFAULT         ; SHERWOOD = 0.037_EB*SC_FILM**ONTH*RE_L**0.8_EB
            CASE(SURF_SPHERICAL) ; SHERWOOD = 2._EB + 0.6_EB*SC_FILM**ONTH*SQRT(RE_L)
         END SELECT        
         
         H_MASS   = MAX(H_MASS_DNS,SHERWOOD*D_FILM/LENGTH_SCALE)
      ENDIF
   ENDIF H_MASS_IF

ENDIF IF_DO_EVAPORATION


! Calculate reaction rates for liquids, solids and vegetation

MATERIAL_LOOP: DO N=1,N_MATS  ! Loop over all materials in the cell (alpha subscript)

   IF (RHO_S(N) < TWO_EPSILON_EB) CYCLE MATERIAL_LOOP  ! If component alpha density is zero, go on to the next material.
   ML => MATERIAL(MATL_INDEX(N))

!   RHO_DOT_REAC = 0._EB
!   RHO_DOT_OUT_REAC = 0._EB
!   M_DOT_G_PPP_ADJUST_REAC = 0._EB
!   M_DOT_G_PPP_ACTUAL_REAC = 0._EB
!   M_DOT_S_PPP_REAC = 0._EB
!   Q_DOT_S_PPP_REAC = 0._EB
!   Q_DOT_G_PPP_REAC = 0._EB
!   Q_DOT_O2_PPP_REAC = 0._EB
!   Q_DOT_PART_REAC = 0._EB
!   M_DOT_PART_REAC = 0._EB
   
   REACTION_LOOP_1: DO J=1,ML%N_REACTIONS  ! Loop over the reactions (beta subscript)

      SELECT CASE (ML%PYROLYSIS_MODEL)

         CASE (PYROLYSIS_LIQUID)

            ! Limit the burning rate to (200 kW/m2) / h_g

             MFLUX_MAX = 200.E3_EB/ML%H_R(J,INT(TMP_F))

            ! Calculate the mass flux of liquid component N at the surface if this is a surface cell.

            IF (DO_EVAPORATION) THEN
               IF (B_NUMBER>TWO_EPSILON_EB) THEN
                  MFLUX = MAX(0._EB,MIN(MFLUX_MAX,H_MASS*RHO_FILM*LOG(1._EB+B_NUMBER)*(Y_SV(N) + (Y_SV(N)-Y_GAS(N))/B_NUMBER)))
               ELSE
                  MFLUX = 0._EB
               ENDIF
            ELSE
               MFLUX = 0._EB
            ENDIF

            IF (DX_S(SOLID_CELL_INDEX)>TWO_EPSILON_EB) THEN

               ! If the liquid temperature (TMP_S) is greater than the boiling temperature of the current liquid component
               ! ((T_BOIL_EFF(N)), calculate the additional mass loss rate of this component (RHO_DOT_EXTRA) necessary to bring
               ! the liquid temperature back to the boiling temperature.

               RHO_DOT_EXTRA = 0._EB
               IF (TMP_S>T_BOIL_EFF(N)) THEN
                  ITMP = MIN(I_MAX_TEMP,INT(TMP_S))
                  H_S = ML%H(ITMP) + (TMP_S-REAL(ITMP,EB))*(ML%H(ITMP+1)-ML%H(ITMP))
                  ITMP = INT(T_BOIL_EFF(N))
                  H_S = H_S - (ML%H(ITMP) + (T_BOIL_EFF(N)-REAL(ITMP,EB))*(ML%H(ITMP+1)-ML%H(ITMP)))
                  H_S = H_S * RHO_S(N)
                  H_R = ML%H_R(1,NINT(T_BOIL_EFF(N)))
                  RHO_DOT_EXTRA = H_S/(H_R*DT_BC)  ! kg/m3/s
               ENDIF

               ! Calculate the mass loss rate per unit volume of this liquid component (RHO_DOT)

               SELECT CASE(SF%GEOMETRY)
                  CASE DEFAULT
                     MFLUX = MIN(MFLUX_MAX,MFLUX + RHO_DOT_EXTRA*DX_S(SOLID_CELL_INDEX))
                     RHO_DOT_REAC(J) =MFLUX/DX_S(SOLID_CELL_INDEX) ! kg/m3/s
                  CASE(SURF_SPHERICAL)
                     NWP = SUM(ONE_D%N_LAYER_CELLS(1:SF%N_LAYERS))
                     R_S_0 = SF%INNER_RADIUS + ONE_D%X(NWP) - ONE_D%X(0)
                     R_S_1 = SF%INNER_RADIUS + ONE_D%X(NWP) - ONE_D%X(1)
                     DR = (R_S_0**3-R_S_1**3)/(3._EB*R_S_0**2)
                     MFLUX = MIN(MFLUX_MAX,MFLUX + RHO_DOT_EXTRA*DR)
                     RHO_DOT_REAC(J) = MFLUX/DR
               END SELECT

            ENDIF

         CASE (PYROLYSIS_SOLID)

            ! Reaction rate in 1/s (Tech Guide: r_alpha_beta)

            REACTION_RATE = ML%A(J)*(RHO_S(N))**ML%N_S(J)*EXP(-ML%E(J)/(R0*TMP_S))

            ! power term

            IF (ABS(ML%N_T(J))>=TWO_EPSILON_EB) REACTION_RATE = REACTION_RATE * TMP_S**ML%N_T(J)

            ! Oxidation reaction?

            IF ( ML%N_O2(J)>0._EB .AND. O2_INDEX>0 ) THEN
               ! Calculate oxygen volume fraction at the surface
               X_O2 = SPECIES(O2_INDEX)%RCON*Y_O2_F/RSUM(IIG,JJG,KKG)
               ! Calculate oxygen concentration inside the material, assuming decay function
               X_O2 = X_O2 * EXP(-DEPTH/(TWO_EPSILON_EB+ML%GAS_DIFFUSION_DEPTH(J)))
               REACTION_RATE = REACTION_RATE * X_O2**ML%N_O2(J)
            ENDIF
            REACTION_RATE = MIN(REACTION_RATE,ML%MAX_REACTION_RATE(J))  ! User-specified limit
            RHO_DOT_REAC(J) = REACTION_RATE  ! Tech Guide: rho_s(0)*r_alpha,beta kg/m3/s

         CASE (PYROLYSIS_SURFACE_OXIDATION)

            ! Reaction rate in kg/m2/s
            REACTION_RATE = ML%A(J)*EXP(-ML%E(J)/(R0*TMP_S))

            ! Estimate surface oxygen concentration from mass transport
            TMP_FILM = (TMP_F+TMP(IIG,JJG,KKG))/2._EB
            ! Get oxygen mass fraction
            ZZ_GET(1:N_TRACKED_SPECIES) = MAX(0._EB,ZZ(IIG,JJG,KKG,1:N_TRACKED_SPECIES))
            CALL GET_MASS_FRACTION(ZZ_GET,O2_INDEX,Y_O2)
            CALL GET_SPECIFIC_HEAT(ZZ_GET,CP_FILM,TMP_FILM)
            ! Mass transfer coefficient
            H_MASS = B1%HEAT_TRANS_COEF/CP_FILM
            ! Mass stoichiometric coefficient for oxygen
            NU_O2_CHAR = ML%NU_GAS_M(O2_INDEX,J)
            Y_O2_S = 0._EB
            IF (H_MASS>0._EB) THEN
               Y_O2_S = (SQRT(4._EB*REACTION_RATE/H_MASS*Y_O2+(REACTION_RATE*NU_O2_CHAR/H_MASS)**2._EB + &
                  2._EB*REACTION_RATE*NU_O2_CHAR/H_MASS+1._EB)-REACTION_RATE*NU_O2_CHAR/H_MASS-1) / &
                  (2._EB*REACTION_RATE/H_MASS)
            ENDIF

            ! Compute LENGTH_SCALE: 1/(surface-to-volume ratio)
            IF (SF%BOUNDARY_FUEL_MODEL) THEN
               LENGTH_SCALE = 1._EB/(SF%SURFACE_VOLUME_RATIO(LAYER_INDEX)*SF%PACKING_RATIO(LAYER_INDEX))
            ELSE
               LENGTH_SCALE = SF%INNER_RADIUS + SUM(ONE_D%LAYER_THICKNESS(1:SF%N_LAYERS))
               SELECT CASE(SF%GEOMETRY)
                  CASE(SURF_SPHERICAL)
                     LENGTH_SCALE = LENGTH_SCALE/3._EB
                  CASE DEFAULT
                     LENGTH_SCALE = LENGTH_SCALE/2._EB
               END SELECT
            ENDIF

            REACTION_RATE = Y_O2_S/LENGTH_SCALE*REACTION_RATE
            REACTION_RATE = MIN(REACTION_RATE,ML%MAX_REACTION_RATE(J))  ! User-specified limit
            RHO_DOT_REAC(J) = MAX(REACTION_RATE,0._EB)  ! Tech Guide: rho_s(0)*r_alpha,beta kg/m3/s

      END SELECT
   ENDDO REACTION_LOOP_1

   RHO_DOT_REAC_SUM = SUM(RHO_DOT_REAC(1:ML%N_REACTIONS))
   IF (RHO_DOT_REAC_SUM < TWO_EPSILON_EB .AND. .NOT. REMOVE_LAYER) CYCLE MATERIAL_LOOP

   IF (REMOVE_LAYER) THEN
      RHO_DOT = RHO_S(N)/DT_BC
      ! If layer is being removed but zero reaction rate, apportion mass loss equally over reactions
      IF (RHO_DOT_REAC_SUM < TWO_EPSILON_EB) RHO_DOT_REAC_SUM = RHO_DOT/ML%N_REACTIONS
   ELSE
      RHO_DOT = MIN(RHO_DOT_REAC_SUM,RHO_S(N)/DT_BC)
   ENDIF

   RHO_DOT_REAC(1:ML%N_REACTIONS) = RHO_DOT_REAC(1:ML%N_REACTIONS)*RHO_DOT/RHO_DOT_REAC_SUM

   ! Optional limiting of mass loss rate based on specified fuel burnout time

   IF (SF%MINIMUM_BURNOUT_TIME<1.E5_EB) RHO_DOT = MIN(RHO_DOT,SF%LAYER_DENSITY(LAYER_INDEX)/SF%MINIMUM_BURNOUT_TIME)

   REACTION_LOOP_2:DO J=1,ML%N_REACTIONS

      RHO_DOT_OUT(N) = RHO_DOT_OUT(N) + RHO_DOT_REAC(J)        

      DO NN=1,ML%N_RESIDUE(J) ! Get residue production (alpha' represents the other materials)
         NNN = FINDLOC(MATL_INDEX,ML%RESIDUE_MATL_INDEX(NN,J),1)
         RHO_DOT_OUT(NNN) = RHO_DOT_OUT(NNN) - ML%NU_RESIDUE(NN,J)*RHO_DOT_REAC(J)
         M_DOT_S_PPP(NNN) = M_DOT_S_PPP(NNN) + ML%NU_RESIDUE(NN,J)*RHO_DOT_REAC(J) ! (m_dot_alpha')'''
      ENDDO

      ! Optional variable heat of reaction

      ITMP = MIN(I_MAX_TEMP,NINT(TMP_S))
      H_R = ML%H_R(J,ITMP)

      ! Calculate various energy and mass source terms

      Q_DOT_S_PPP    = Q_DOT_S_PPP    - RHO_DOT_REAC(J)*H_R ! Tech Guide: q_dot_s'''
      M_DOT_S_PPP(N) = M_DOT_S_PPP(N) - RHO_DOT_REAC(J)     ! m_dot_alpha''' = -rho_s(0) * sum_beta r_alpha,beta
      TMP_G = TMP(IIG,JJG,KKG)
      ! Tech Guide: m_dot_gamma'''

      M_DOT_G_PPP_ACTUAL(:) = M_DOT_G_PPP_ACTUAL(:) + ML%NU_GAS(:,J)*RHO_DOT_REAC(J)
      M_DOT_G_PPP_ADJUST(:) = M_DOT_G_PPP_ADJUST(:) + ML%NU_GAS(:,J)*RHO_DOT_REAC(J) * ML%ADJUST_BURN_RATE(:,J)
      DO NS=1,N_SPECIES
         IF (ML%NU_GAS_P(NS,J) <= 0._EB .AND. ML%NU_GAS_M(NS,J) <= 0._EB) CYCLE
         CALL INTERPOLATE1D_UNIFORM(0,SPECIES(NS)%H_G,TMP_S,H_S_B)
         CALL INTERPOLATE1D_UNIFORM(0,SPECIES(NS)%H_G,TMP_G,H_S)
         IF (ML%NU_GAS_P(NS,J) > 0._EB) &
            Q_DOT_G_PPP = Q_DOT_G_PPP + ML%ADJUST_BURN_RATE_P(NS,J)*ML%NU_GAS_P(NS,J)*RHO_DOT_REAC(J)*(H_S-H_S_B)
         IF (ML%NU_GAS_M(NS,J) > 0._EB) &
            Q_DOT_G_PPP = Q_DOT_G_PPP -                             ML%NU_GAS_M(NS,J)*RHO_DOT_REAC(J)*(H_S-H_S_B)
      ENDDO

      IF (ANY(ML%NU_LPC(:,J)>0._EB)) THEN
         DO NP=1,ML%N_LPC(J)
            IF (ML%NU_LPC(NP,J)<=0._EB) CYCLE
            NP2 = ML%LPC_INDEX(NP,J)
            IF (SF%MATL_PART_INDEX(NP2)==NP) THEN
               M_DOT_PART(NP2)=ML%NU_LPC(NP,J)*RHO_DOT_REAC(J)
               Q_DOT_PART(NP2)=GET_PARTICLE_ENTHALPY(NP,TMP_S)*M_DOT_PART(NP2)
            ENDIF
         ENDDO
      ENDIF

      ! If there is char oxidation, save the HRR per unit volume generated

      IF (ML%N_O2(J)>0._EB) Q_DOT_O2_PPP = Q_DOT_O2_PPP - RHO_DOT_REAC(J)*H_R

   ENDDO REACTION_LOOP_2
   
ENDDO MATERIAL_LOOP

END SUBROUTINE PYROLYSIS


!> \brief Compute the convective heat transfer coefficient

REAL(EB) FUNCTION HEAT_TRANSFER_COEFFICIENT(NMX,DELTA_N_TMP,H_FIXED,SFX,WALL_INDEX_IN,CFACE_INDEX_IN,PARTICLE_INDEX_IN)

USE TURBULENCE, ONLY: LOGLAW_HEAT_FLUX_MODEL,NATURAL_CONVECTION_MODEL,FORCED_CONVECTION_MODEL,RAYLEIGH_HEAT_FLUX_MODEL,&
                      FM_HEAT_FLUX_MODEL
USE PHYSICAL_FUNCTIONS, ONLY: GET_CONDUCTIVITY,GET_VISCOSITY,GET_SPECIFIC_HEAT
REAL(EB), INTENT(IN) :: DELTA_N_TMP,H_FIXED
INTEGER, INTENT(IN) :: NMX
INTEGER, INTENT(IN), OPTIONAL :: WALL_INDEX_IN,PARTICLE_INDEX_IN,CFACE_INDEX_IN
INTEGER  :: SURF_GEOMETRY,ITMP,I,HTR
REAL(EB) :: RE,H_NATURAL,H_FORCED,FRICTION_VELOCITY=0._EB,YPLUS=0._EB,ZSTAR,DN,TMP_FILM,MU_G,K_G,CP_G,&
            TMP_G,R_DROP,RHO_G,ZZ_G(1:N_TRACKED_SPECIES),CONV_LENGTH,GR,RA,NUSSELT_FORCED,NUSSELT_FREE,NUSSELT_IMPINGE,&
            FILM_FAC,PHI,C0_IMP,C1_IMP,C2_IMP,M_IMP
INTEGER, PARAMETER :: NATURAL=1,FORCED=2,IMPACT=3,RESOLVED=4
TYPE(MESH_TYPE), POINTER :: MX
TYPE(SURFACE_TYPE), INTENT(IN), POINTER :: SFX
TYPE(WALL_TYPE), POINTER :: WCX
TYPE(CFACE_TYPE), POINTER :: CFAX
TYPE(LAGRANGIAN_PARTICLE_TYPE), POINTER :: LPX
TYPE(BOUNDARY_PROP1_TYPE), POINTER :: P1X
TYPE(BOUNDARY_PROP2_TYPE), POINTER :: P2X
TYPE(BOUNDARY_COORD_TYPE), POINTER :: BCX

MX  => MESHES(NMX)
CONV_LENGTH = SFX%CONV_LENGTH

! Determine if this is a particle or wall cell

IF (PRESENT(PARTICLE_INDEX_IN)) THEN
   LPX => LAGRANGIAN_PARTICLE(PARTICLE_INDEX_IN)
   P1X => BOUNDARY_PROP1(LPX%B1_INDEX)
   BCX => BOUNDARY_COORD(LPX%BC_INDEX)
   DN = SFX%CONV_LENGTH
   R_DROP = LPX%RADIUS
   IF (R_DROP>TWO_EPSILON_EB) CONV_LENGTH = 2._EB*R_DROP
   TMP_G  = P1X%TMP_G
   RHO_G  = P1X%RHO_G
   ZZ_G(1:N_TRACKED_SPECIES) = P1X%ZZ_G(1:N_TRACKED_SPECIES)
ELSEIF (PRESENT(WALL_INDEX_IN)) THEN
   WCX   => MX%WALL(WALL_INDEX_IN)
   P1X => MX%BOUNDARY_PROP1(WCX%B1_INDEX)
   P2X => MX%BOUNDARY_PROP2(WCX%B2_INDEX)
   BCX => MX%BOUNDARY_COORD(WCX%BC_INDEX)
   DN = 1._EB/P1X%RDN
   TMP_G  = P1X%TMP_G
   RHO_G  = P1X%RHO_G
   ZZ_G(1:N_TRACKED_SPECIES) = P1X%ZZ_G(1:N_TRACKED_SPECIES)
ELSEIF (PRESENT(CFACE_INDEX_IN)) THEN
   CFAX => CFACE(CFACE_INDEX_IN)
   P1X => BOUNDARY_PROP1(CFAX%B1_INDEX)
   P2X => BOUNDARY_PROP2(CFAX%B2_INDEX)
   BCX => BOUNDARY_COORD(CFAX%BC_INDEX)
   DN = 1._EB/P1X%RDN
   TMP_G  = P1X%TMP_G
   RHO_G  = P1X%RHO_G
   ZZ_G(1:N_TRACKED_SPECIES) = P1X%ZZ_G(1:N_TRACKED_SPECIES)
ELSE
   IF (H_FIXED >= 0._EB) THEN
      HEAT_TRANSFER_COEFFICIENT = H_FIXED
   ELSE
      HEAT_TRANSFER_COEFFICIENT = 1.31_EB*ABS(DELTA_N_TMP)**ONTH  ! Natural convection for vertical plane
   ENDIF
   RETURN
ENDIF

! Calculate HEAT_TRANSFER_COEFFICIENT

SELECT CASE(SFX%GEOMETRY)
   CASE DEFAULT                                  ; FILM_FAC = PLATE_FILM_FAC
   CASE(SURF_SPHERICAL)                          ; FILM_FAC = SPHERE_FILM_FAC
   CASE(SURF_CYLINDRICAL,SURF_INNER_CYLINDRICAL) ; FILM_FAC = PLATE_FILM_FAC
   CASE(SURF_CARTESIAN)                          ; FILM_FAC = PLATE_FILM_FAC
END SELECT

TMP_FILM = P1X%TMP_F+FILM_FAC*(TMP_G-P1X%TMP_F)

! If the user wants a specified HTC, set it and return

H_FIXED_IF: IF (H_FIXED >= 0._EB) THEN

   HEAT_TRANSFER_COEFFICIENT = H_FIXED

ELSE H_FIXED_IF

   ! Either do a DNS /SOLID_PHASE_ONLY calcuation or an LES / non-SOLID_PHASE_ONLY calculation.
   ! If this is a DNS calculation at a solid wall, set HTC and return. If SOLID_PHASE do not return until after blowing check.

   HTC_TYPE_IF: IF ( (SIM_MODE==DNS_MODE .OR. SOLID_PHASE_ONLY) .AND. (PRESENT(WALL_INDEX_IN) .OR. PRESENT(CFACE_INDEX_IN)) ) THEN

      HEAT_TRANSFER_COEFFICIENT = 2._EB * P1X%K_G * P1X%RDN
      IF (SIM_MODE==DNS_MODE) RETURN

   ELSE HTC_TYPE_IF

      CALL GET_VISCOSITY(ZZ_G,MU_G,TMP_FILM)
      CALL GET_CONDUCTIVITY(ZZ_G,K_G,TMP_FILM)

      HTC_MODEL_SELECT: SELECT CASE(SFX%HEAT_TRANSFER_MODEL)
         CASE(DEFAULT_HTC_MODEL,IMPINGING_JET_HTC_MODEL)
            RE = RHO_G*P1X%U_TANG*CONV_LENGTH/MU_G
            GR = GRAV*ABS(DELTA_N_TMP)*CONV_LENGTH**3*(RHO_G/MU_G)**2/TMP_FILM
            IF (SFX%BOUNDARY_FUEL_MODEL) THEN
               SURF_GEOMETRY = SURF_CYLINDRICAL
            ELSE
               SURF_GEOMETRY = SFX%GEOMETRY
            ENDIF
            ! Check if custom Nusselt correlation is defined
            IF (ANY((/SFX%NUSSELT_C0,SFX%NUSSELT_C1,SFX%NUSSELT_C2,SFX%NUSSELT_M/)>0._EB)) THEN
               CALL FORCED_CONVECTION_MODEL(NUSSELT_FORCED,RE,PR_ONTH,SURF_GEOMETRY,&
                  SFX%NUSSELT_C0,SFX%NUSSELT_C1,SFX%NUSSELT_C2,SFX%NUSSELT_M)
            ELSE
               CALL FORCED_CONVECTION_MODEL(NUSSELT_FORCED,RE,PR_ONTH,SURF_GEOMETRY)
            ENDIF
            RA = GR*PR_AIR
            CALL NATURAL_CONVECTION_MODEL(NUSSELT_FREE,RA,SFX,BCX%IOR,DELTA_N_TMP)
            NUSSELT_IMPINGE = 0._EB
            IF (SFX%HEAT_TRANSFER_MODEL==IMPINGING_JET_HTC_MODEL) THEN
               ! Huang, G.C. : Investigations of Heat-Transfer Coefficients for Air Flow Through Round Jets Impinging
               ! Normal to a Heat-Transfer Surface. J. Heat Transfer, vol. 85, no. 3, Aug. 1963, pp. 237-245.
               RE = RHO_G*P1X%U_IMPACT*CONV_LENGTH/MU_G
               C0_IMP = 0.000_EB; IF (SFX%NUSSELT_C0>0._EB) THEN; C0_IMP = SFX%NUSSELT_C0; ENDIF
               C1_IMP = 0.055_EB; IF (SFX%NUSSELT_C1>0._EB) THEN; C1_IMP = SFX%NUSSELT_C1; ENDIF
               C2_IMP = 0.000_EB; IF (SFX%NUSSELT_C2>0._EB) THEN; C2_IMP = SFX%NUSSELT_C2; ENDIF
               M_IMP  = 0.800_EB; IF (SFX%NUSSELT_M >0._EB) THEN; M_IMP  = SFX%NUSSELT_M ; ENDIF
               NUSSELT_IMPINGE = C0_IMP + (C1_IMP*RE**M_IMP - C2_IMP) * PR_ONTH
            ENDIF
            IF (PRESENT(PARTICLE_INDEX_IN)) THEN
               HEAT_TRANSFER_COEFFICIENT = MAX(NUSSELT_FORCED,NUSSELT_FREE)*K_G/CONV_LENGTH
               ! BOUNDARY_PROP2 not allocated for particles
            ELSE
               P2X%HEAT_TRANSFER_REGIME=MAXLOC((/NUSSELT_FREE,NUSSELT_FORCED,NUSSELT_IMPINGE,2._EB*CONV_LENGTH/DN/),DIM=1)
               SELECT CASE(P2X%HEAT_TRANSFER_REGIME)
                  CASE(NATURAL)
                     P2X%Z_STAR = NUSSELT_FREE**0.25_EB
                  CASE(FORCED)
                     P2X%Z_STAR = SQRT(NUSSELT_FORCED)
                  CASE(IMPACT)
                     P2X%Z_STAR = SQRT(NUSSELT_IMPINGE)
                  CASE(RESOLVED)
                     P2X%Z_STAR = 1._EB
               END SELECT
               HEAT_TRANSFER_COEFFICIENT = MAX(NUSSELT_FREE,NUSSELT_FORCED,NUSSELT_IMPINGE,2._EB*CONV_LENGTH/DN)*K_G/CONV_LENGTH
            ENDIF
         CASE(LOGLAW_HTC_MODEL)
            CALL GET_SPECIFIC_HEAT(ZZ_G,CP_G,TMP_FILM)
            FRICTION_VELOCITY = P2X%U_TAU
            YPLUS = P2X%Y_PLUS
            CALL LOGLAW_HEAT_FLUX_MODEL(H_FORCED,YPLUS,FRICTION_VELOCITY,K_G,RHO_G,CP_G,MU_G)
            HEAT_TRANSFER_COEFFICIENT = H_FORCED
         CASE(RAYLEIGH_HTC_MODEL)
            ! not appropriate for a particle, used with SURF and CFACE only
            CALL GET_SPECIFIC_HEAT(ZZ_G,CP_G,TMP_FILM)
            CALL RAYLEIGH_HEAT_FLUX_MODEL(H_NATURAL,ZSTAR,HTR,DN,P1X%TMP_F,TMP_G,K_G,RHO_G,CP_G,MU_G,P1X%U_TANG,P1X%U_IMPACT)
            P2X%Z_STAR = ZSTAR
            P2X%HEAT_TRANSFER_REGIME = HTR
            HEAT_TRANSFER_COEFFICIENT = H_NATURAL
         CASE(FM_HTC_MODEL)
            CALL GET_SPECIFIC_HEAT(ZZ_G,CP_G,TMP_FILM)
            CALL FM_HEAT_FLUX_MODEL(H_NATURAL,DN,P1X%TMP_F,TMP_G,K_G,RHO_G,CP_G)
            HEAT_TRANSFER_COEFFICIENT = H_NATURAL
      END SELECT HTC_MODEL_SELECT
   ENDIF HTC_TYPE_IF
ENDIF H_FIXED_IF

IF (SFX%BLOWING .AND. .NOT. SFX%BOUNDARY_FUEL_MODEL .AND. SIM_MODE /= DNS_MODE .AND. ALLOCATED(P1X%M_DOT_G_PP_ACTUAL)) THEN
   PHI = 0._EB
   ITMP = INT(TMP_FILM)
   DO I=1,N_TRACKED_SPECIES
      IF (ABS(P1X%M_DOT_G_PP_ACTUAL(I)) <= TWO_EPSILON_EB) CYCLE
      PHI = PHI + P1X%M_DOT_G_PP_ACTUAL(I)*CP_Z(ITMP,I)
   ENDDO
   IF (ABS(PHI)>TWO_EPSILON_EB .AND. ABS(HEAT_TRANSFER_COEFFICIENT)>=TWO_EPSILON_EB) THEN
      PHI = PHI / HEAT_TRANSFER_COEFFICIENT
      IF (PHI > 100._EB) THEN
         HEAT_TRANSFER_COEFFICIENT = 0._EB
      ELSE
         HEAT_TRANSFER_COEFFICIENT = HEAT_TRANSFER_COEFFICIENT * PHI/(EXP(PHI)-1._EB)
      ENDIF
   ENDIF
ENDIF

END FUNCTION HEAT_TRANSFER_COEFFICIENT


!> \brief Transfer heat from one 3-D sweep dirction to the other two.
!> \param NM Mesh index
!> \details Sweep over all 3D wall cells and look for those that are not oriented in the current sweep direction.
!> Copy the updated temperatures of the sweep direction to the non-sweep directions.

SUBROUTINE HT3D_TEMPERATURE_EXCHANGE(NM)

USE COMP_FUNCTIONS, ONLY: CURRENT_TIME
INTEGER, INTENT(IN) :: NM
INTEGER :: NWP,I_NODE,I,II,IWA,NM2,IW,ITW
REAL(EB) :: TMP_1,TMP_NWP,TNOW
TYPE(BOUNDARY_ONE_D_TYPE), POINTER :: ONE_D2
TYPE(BOUNDARY_THR_D_TYPE), POINTER :: THR_D
TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC,BC2
TYPE(WALL_TYPE), POINTER :: WC,WC2
TYPE(THIN_WALL_TYPE), POINTER :: TW,TW2
TYPE(MESH_TYPE), POINTER :: M
TYPE(SURFACE_TYPE), POINTER :: SF
TYPE(BOUNDARY_ONE_D_TYPE), POINTER :: ONE_D

TNOW = CURRENT_TIME()

M => MESHES(NM)

!$OMP PARALLEL

!$OMP DO SCHEDULE(GUIDED) PRIVATE(IW,WC,SF,BC,ONE_D,NWP,THR_D,TMP_1,TMP_NWP,I,II,IWA,NM2,I_NODE,WC2,BC2,ONE_D2,TW2)
WALL_LOOP: DO IW=1,M%N_EXTERNAL_WALL_CELLS+M%N_INTERNAL_WALL_CELLS

   WC => M%WALL(IW)
   SF => SURFACE(WC%SURF_INDEX)
   IF (SF%HT_DIM==1 .OR. WC%BOUNDARY_TYPE/=SOLID_BOUNDARY) CYCLE WALL_LOOP
   BC => M%BOUNDARY_COORD(WC%BC_INDEX)
   IF (ABS(BC%IOR)==M%HT_3D_SWEEP_DIRECTION) CYCLE WALL_LOOP

   ONE_D => M%BOUNDARY_ONE_D(WC%OD_INDEX)
   NWP = SUM(ONE_D%N_LAYER_CELLS(1:ONE_D%N_LAYERS))
   THR_D => M%BOUNDARY_THR_D(WC%TD_INDEX)

   TMP_1   = ONE_D%TMP(1)
   TMP_NWP = ONE_D%TMP(NWP)
   NODE_LOOP: DO I=1,NWP  ! Nodes of the chosen wall cell.
      IF (.NOT.THR_D%NODE(I)%HT3D) CYCLE NODE_LOOP
      IF (THR_D%NODE(I)%ALTERNATE_WALL_COUNT==0) CYCLE NODE_LOOP
      IF (.NOT.ANY(ABS(THR_D%NODE(I)%ALTERNATE_WALL_IOR(:))==M%HT_3D_SWEEP_DIRECTION)) CYCLE NODE_LOOP
      WEIGHT_LOOP: DO II=1,THR_D%NODE(I)%ALTERNATE_WALL_COUNT
         IWA = THR_D%NODE(I)%ALTERNATE_WALL_INDEX(II)
         NM2 = THR_D%NODE(I)%ALTERNATE_WALL_MESH(II)
         I_NODE = THR_D%NODE(I)%ALTERNATE_WALL_NODE(II)
         IF (THR_D%NODE(I)%ALTERNATE_WALL_TYPE(II)==0) THEN
            WC2 => MESHES(NM2)%WALL(IWA)
            BC2 => MESHES(NM2)%BOUNDARY_COORD(WC2%BC_INDEX)
            ONE_D2 => MESHES(NM2)%BOUNDARY_ONE_D(WC2%OD_INDEX)
         ELSE
            TW2 => MESHES(NM2)%THIN_WALL(IWA)
            BC2 => MESHES(NM2)%BOUNDARY_COORD(TW2%BC_INDEX)
            ONE_D2 => MESHES(NM2)%BOUNDARY_ONE_D(TW2%OD_INDEX)
         ENDIF
         IF (ABS(BC2%IOR)/=M%HT_3D_SWEEP_DIRECTION) CYCLE WEIGHT_LOOP
         ONE_D%TMP(I) = (ONE_D%RHO_C_S(I)*ONE_D%TMP(I) + &
                         THR_D%NODE(I)%ALTERNATE_WALL_WEIGHT(II)*ONE_D2%RHO_C_S(I_NODE)*ONE_D2%DELTA_TMP(I_NODE))/ONE_D%RHO_C_S(I)
      ENDDO WEIGHT_LOOP
      IF (THR_D%NODE(I)%MESH_NUMBER==NM) M%TMP(THR_D%NODE(I)%I,THR_D%NODE(I)%J,THR_D%NODE(I)%K) = ONE_D%TMP(I)
   ENDDO NODE_LOOP
   ONE_D%TMP(0)     = ONE_D%TMP(0)     + ONE_D%TMP(1)   - TMP_1
   ONE_D%TMP(NWP+1) = ONE_D%TMP(NWP+1) + ONE_D%TMP(NWP) - TMP_NWP

ENDDO WALL_LOOP
!$OMP END DO

!$OMP DO SCHEDULE(GUIDED) PRIVATE(ITW,TW,SF,BC,ONE_D,NWP,THR_D,I,II,IWA,NM2,I_NODE,WC2,BC2,ONE_D2,TW2)
THIN_WALL_LOOP: DO ITW=1,M%N_THIN_WALL_CELLS

   TW => M%THIN_WALL(ITW)
   SF => SURFACE(TW%SURF_INDEX)
   BC => M%BOUNDARY_COORD(TW%BC_INDEX)
   IF (ABS(BC%IOR)==M%HT_3D_SWEEP_DIRECTION) CYCLE THIN_WALL_LOOP
   ONE_D => M%BOUNDARY_ONE_D(TW%OD_INDEX)
   NWP = SUM(ONE_D%N_LAYER_CELLS(1:SF%N_LAYERS))
   THR_D => M%BOUNDARY_THR_D(TW%TD_INDEX)
   IF (.NOT.ALLOCATED(THR_D%NODE)) CYCLE THIN_WALL_LOOP

   TMP_1   = ONE_D%TMP(1)
   TMP_NWP = ONE_D%TMP(NWP)
   NODE_LOOP_2: DO I=1,NWP  ! Nodes of the chosen wall cell.
      IF (THR_D%NODE(I)%ALTERNATE_WALL_COUNT==0) CYCLE NODE_LOOP_2
      IF (.NOT.ANY(ABS(THR_D%NODE(I)%ALTERNATE_WALL_IOR(:))==M%HT_3D_SWEEP_DIRECTION)) CYCLE NODE_LOOP_2
      WEIGHT_LOOP_2: DO II=1,THR_D%NODE(I)%ALTERNATE_WALL_COUNT
         IWA = THR_D%NODE(I)%ALTERNATE_WALL_INDEX(II)
         NM2 = THR_D%NODE(I)%ALTERNATE_WALL_MESH(II)
         I_NODE = THR_D%NODE(I)%ALTERNATE_WALL_NODE(II)
         IF (THR_D%NODE(I)%ALTERNATE_WALL_TYPE(II)==0) THEN
            WC2 => MESHES(NM2)%WALL(IWA)
            BC2 => MESHES(NM2)%BOUNDARY_COORD(WC2%BC_INDEX)
            ONE_D2 => MESHES(NM2)%BOUNDARY_ONE_D(WC2%OD_INDEX)
         ELSE
            TW2 => MESHES(NM2)%THIN_WALL(IWA)
            BC2 => MESHES(NM2)%BOUNDARY_COORD(TW2%BC_INDEX)
            ONE_D2 => MESHES(NM2)%BOUNDARY_ONE_D(TW2%OD_INDEX)
         ENDIF
         IF (ABS(BC2%IOR)/=M%HT_3D_SWEEP_DIRECTION) CYCLE WEIGHT_LOOP_2
         ONE_D%TMP(I) = ONE_D%TMP(I) + THR_D%NODE(I)%ALTERNATE_WALL_WEIGHT(II)*ONE_D2%DELTA_TMP(I_NODE)
      ENDDO WEIGHT_LOOP_2
   ENDDO NODE_LOOP_2
   ONE_D%TMP(0)     = ONE_D%TMP(0)     + ONE_D%TMP(1)   - TMP_1
   ONE_D%TMP(NWP+1) = ONE_D%TMP(NWP+1) + ONE_D%TMP(NWP) - TMP_NWP

ENDDO THIN_WALL_LOOP
!$OMP END DO

!$OMP END PARALLEL

T_USED(6)=T_USED(6)+CURRENT_TIME()-TNOW
END SUBROUTINE HT3D_TEMPERATURE_EXCHANGE


!> \brief Perform a numerical TGA (thermo-gravimetric analysis) at the start of the simulation
!> \param NM Mesh index

SUBROUTINE TGA_ANALYSIS(NM)

USE PHYSICAL_FUNCTIONS, ONLY: SURFACE_DENSITY
USE COMP_FUNCTIONS, ONLY: SHUTDOWN
REAL(EB) :: DT_TGA=0.01_EB,T_TGA,SURF_DEN_0,HRR
REAL(EB), ALLOCATABLE, DIMENSION(:) :: SURF_DEN
INTEGER, INTENT(IN) :: NM
INTEGER :: N_TGA,I,IW,IP,N,NR
CHARACTER(80) :: MESSAGE,TCFORM
TYPE(SURFACE_TYPE), POINTER :: SF
TYPE(WALL_TYPE), POINTER :: WC
TYPE(LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP
TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1
TYPE(BOUNDARY_ONE_D_TYPE), POINTER :: ONE_D

CALL POINT_TO_MESH(NM)

SF => SURFACE(TGA_SURF_INDEX)
ALLOCATE(SURF_DEN(0:SF%N_MATL))
RADIATION = .FALSE.
TGA_HEATING_RATE = TGA_HEATING_RATE/60._EB  ! K/min --> K/s
TGA_FINAL_TEMPERATURE = TGA_FINAL_TEMPERATURE + TMPM  ! C --> K
N_TGA = NINT((TGA_FINAL_TEMPERATURE-TMPA)/(TGA_HEATING_RATE*DT_TGA))
T_TGA = 0._EB

IF (TGA_WALL_INDEX>0) THEN
   IW = TGA_WALL_INDEX
   WC => WALL(IW)
   ONE_D => BOUNDARY_ONE_D(WC%OD_INDEX)
   B1 => BOUNDARY_PROP1(WC%B1_INDEX)
ELSEIF (TGA_PARTICLE_INDEX>0) THEN
   IP = TGA_PARTICLE_INDEX
   LP => LAGRANGIAN_PARTICLE(IP)
   ONE_D => BOUNDARY_ONE_D(LP%OD_INDEX)
   B1 => BOUNDARY_PROP1(LP%B1_INDEX)
ELSE
   WRITE(MESSAGE,'(3A)') 'ERROR(370): SURF ',TRIM(SF%ID),' No wall or particle for TGA_ANALYSIS.'
   CALL SHUTDOWN(MESSAGE) ; RETURN
ENDIF

OPEN (LU_TGA,FILE=FN_TGA,FORM='FORMATTED',STATUS='REPLACE')
WRITE(TCFORM,'(A,I3.3,A,I3.3,A)') "(A,",SF%N_MATL+1,"(A,',')",SF%N_MATL+1,"(A,','),A)"
WRITE(LU_TGA,TCFORM) 's,C,',('g/g',N=1,SF%N_MATL+1),('1/s',N=1,SF%N_MATL+1),'W/g,W/g'
WRITE(LU_TGA,TCFORM) 'Time,Temp,','Total Mass',(TRIM(MATERIAL(ONE_D%MATL_INDEX(N))%ID)//' Mass',N=1,SF%N_MATL),'Total MLR',&
                     (TRIM(MATERIAL(ONE_D%MATL_INDEX(N))%ID)//' MLR',N=1,SF%N_MATL),'MCC,DSC'

SURF_DEN_0 = SF%SURFACE_DENSITY
WRITE(TCFORM,'(A,I3.3,5A)') "(",2*SF%N_MATL+5,"(",TRIM(FMT_R),",','),",TRIM(FMT_R),")"

DO I=1,N_TGA
   IF (ONE_D%LAYER_THICKNESS(1)<TWO_EPSILON_EB) EXIT
   T_TGA = I*DT_TGA
   B1%TMP_G = TMPA + TGA_HEATING_RATE*T_TGA
   IF (TGA_WALL_INDEX>0) THEN
      CALL SOLID_HEAT_TRANSFER(NM,T_TGA,DT_TGA,WALL_INDEX=IW)
   ELSE
      CALL SOLID_HEAT_TRANSFER(NM,T_TGA,DT_TGA,PARTICLE_INDEX=IP)
   ENDIF
   IF (I==1 .OR. MOD(I,NINT(1._EB/(TGA_HEATING_RATE*DT_TGA)))==0) THEN
      IF (TGA_WALL_INDEX>0) THEN
         SURF_DEN(0) = SURFACE_DENSITY(0,SF,ONE_D)
         DO N=1,SF%N_MATL
            SURF_DEN(N) = SURFACE_DENSITY(0,SF,ONE_D,MATL_INDEX=N)
         ENDDO
      ELSE
         SURF_DEN(0) = SURFACE_DENSITY(0,SF,ONE_D)
         DO N=1,SF%N_MATL
            SURF_DEN(N) = SURFACE_DENSITY(0,SF,ONE_D,MATL_INDEX=N)
         ENDDO
      ENDIF
      IF (N_REACTIONS>0) THEN
         HRR = 0._EB
         DO NR=1,N_REACTIONS
            HRR = B1%M_DOT_G_PP_ADJUST(REACTION(NR)%FUEL_SMIX_INDEX)*0.001*REACTION(NR)%HEAT_OF_COMBUSTION / &
                  (B1%AREA_ADJUST*SURF_DEN_0)
         ENDDO
      ELSE
         HRR = 0._EB
      ENDIF
      WRITE(LU_TGA,TCFORM) REAL(T_TGA,FB), REAL(B1%TMP_G-TMPM,FB), (REAL(SURF_DEN(N)/SURF_DEN_0,FB),N=0,SF%N_MATL), &
                           REAL(-SUM(ONE_D%M_DOT_S_PP(1:SF%N_MATL))/SURF_DEN_0,FB), &
                           (REAL(-ONE_D%M_DOT_S_PP(N)/SURF_DEN_0,FB),N=1,SF%N_MATL), &
                           REAL(HRR,FB), REAL(B1%HEAT_TRANS_COEF*(B1%TMP_G-B1%TMP_F)*0.001_EB/SURF_DEN_0,FB)
   ENDIF
ENDDO

CLOSE(LU_TGA)
DEALLOCATE(SURF_DEN)

END SUBROUTINE TGA_ANALYSIS

END MODULE WALL_ROUTINES

